{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tutorial: VIC 5 Image Driver Parameter Conversion\n",
    "====\n",
    "***Converting parameters from ASCII VIC 4 format to netCDF VIC 5 Image Driver Format***\n",
    "\n",
    "This Jupyter Notebook outlines one approach to converting VIC parameters from ASCII to netCDF format. For this tutorial, we'll convert three datasets from ASCII to netCDF:\n",
    "\n",
    "1. **Livneh et al. (2015) - 1/16th deg. VIC parameters**\n",
    "  - Description: http://www.colorado.edu/lab/livneh/data\n",
    "  - Data: ftp://livnehpublicstorage.colorado.edu/public/Livneh.2015.NAmer.Dataset/nldas.vic.params/\n",
    "  - Citation: Livneh B., E.A. Rosenberg, C. Lin, B. Nijssen, V. Mishra, K.M. Andreadis, E.P. Maurer, and D.P. Lettenmaier, 2013: A Long-term hydrologically based dataset of land surface fluxes and states for the conterminous United States: update and extensions, Journal of Climate, 26, 9384–9392.\n",
    "2. **Global 1/2 deg. VIC parameters**\n",
    "  - Description: http://www.hydro.washington.edu/SurfaceWaterGroup/Data/vic_global_0.5deg.html\n",
    "  - Data: ftp://ftp.hydro.washington.edu/pub/HYDRO/data/VIC_param/vic_params_global_0.5deg.tgz\n",
    "  - Citation: Nijssen, B.N., G.M. O'Donnell, D.P. Lettenmaier and E.F. Wood, 2001: Predicting the discharge of global rivers, J. Clim., 14(15), 3307-3323, doi: 10.1175/1520-0442(2001)014<3307:PTDOGR>2.0.CO;2.\n",
    "3. **Maurer et al. (2002) - 1/8th deg. VIC parameters**\n",
    "  - Description: http://www.hydro.washington.edu/SurfaceWaterGroup/Data/VIC_retrospective/index.html\n",
    "  - Data: http://www.hydro.washington.edu/SurfaceWaterGroup/Data/VIC_retrospective/index.html\n",
    "  - Citation: Maurer, E.P., A.W. Wood, J.C. Adam, D.P. Lettenmaier, and B. Nijssen, 2002: A long-term hydrologically-based data set of land surface fluxes and states for the conterminous United States, J. Climate 15, 3237-3251. \n",
    "\n",
    "All of these datasets include the following parameter sets:\n",
    "- Soil Parameter file\n",
    "- Vegetation Library file\n",
    "- Vegetation Parameter file\n",
    "- Snowbands file \n",
    "\n",
    "### Outputs\n",
    "For each of the parameter sets above, we'll be producing two files:\n",
    "1.  VIC 5 Image Driver Input Parameters (netCDF file defining model parameters)\n",
    "2.  VIC 5 Image Driver Domain File (netCDF file defining spatial extent of model domain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python Imports and Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "python version : 3.5.1 |Continuum Analytics, Inc.| (default, Dec  7 2015, 11:24:55) \n",
      "[GCC 4.2.1 (Apple Inc. build 5577)]\n",
      "numpy version  : 1.10.4\n",
      "xarray version : 0.7.2\n",
      "User           : jhamman\n",
      "Date/Time      : 2016-06-08 15:01:23.814110\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import getpass\n",
    "from datetime import datetime\n",
    "\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# For more information on tonic, see: https://github.com/UW-Hydro/tonic/\n",
    "import tonic.models.vic.grid_params as gp\n",
    "\n",
    "# Metadata to be used later\n",
    "user = getpass.getuser()\n",
    "now = datetime.now()\n",
    "\n",
    "print('python version : %s' % os.sys.version)\n",
    "print('numpy version  : %s' % np.version.full_version)\n",
    "print('xarray version : %s' % xr.version.version)\n",
    "print('User           : %s' % user)\n",
    "print('Date/Time      : %s' % now)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set Path Information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set the path to the datasets here\n",
    "dpath = './'  # root input data path\n",
    "opath = './'  # output data path\n",
    "ldpath = os.path.join(dpath, 'Livneh_0.0625_NLDAS')  # Path to Livneh Parameters\n",
    "gdpath = os.path.join(dpath, 'Nijssen_0.5_Global')  # Path to Global Parameters\n",
    "mdpath = os.path.join(dpath, 'Maurer_0.125_NLDAS')  # Path to Maurer Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Livneh Domain File\n",
    "\n",
    "Along with the VIC model parameters, we also need a domain file that describes the spatial extent and active grid cells in the model domain. The domain file must exactly match the parameters and forcings. The Livneh dataset includes a DEM file in netCDF format that we can use to construct the domain file.\n",
    "\n",
    "Steps:\n",
    "\n",
    "1. Open dem\n",
    "2. Set useful global attributes\n",
    "3. Create the mask/frac variables using the non-missing dem mask.\n",
    "4. Calculate the grid cell area using `cdo`\n",
    "5. Add the grid cell area back into the domain dataset.\n",
    "6. Save the domain dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dom_file = os.path.join(opath, 'domain.vic.global0.0625deg.%s.nc' % now.strftime('%Y%m%d'))\n",
    "dem = xr.open_dataset(os.path.join(ldpath, 'Composite.DEM.NLDAS.mex.0625.nc'))\n",
    "\n",
    "dom_ds = xr.Dataset()\n",
    "\n",
    "# Set global attributes\n",
    "dom_ds.attrs['title'] = 'VIC domain data'\n",
    "dom_ds.attrs['Conventions'] = 'CF-1.6'\n",
    "dom_ds.attrs['history'] = 'created by %s, %s' % (user, now)\n",
    "dom_ds.attrs['user_comment'] = 'VIC domain data'\n",
    "dom_ds.attrs['source'] = 'generated from VIC North American 1/16 deg. model parameters, see Livneh et al. (2015) for more information'\n",
    "\n",
    "# since we have it, put the elevation in the domain file\n",
    "dom_ds['elev'] = dem['Band1']\n",
    "dom_ds['elev'].attrs['long_name'] = 'gridcell_elevation'\n",
    "dom_ds['elev'].attrs['units'] = 'm'\n",
    "\n",
    "# Get the mask variable\n",
    "dom_ds['mask'] = dem['Band1'].notnull().astype(np.int)\n",
    "dom_ds['mask'].attrs['long_name'] = 'domain mask'\n",
    "dom_ds['mask'].attrs['comment'] = '0 indicates cell is not active'\n",
    "\n",
    "# For now, the frac variable is going to be just like the mask\n",
    "dom_ds['frac'] = dom_ds['mask'].astype(np.float)\n",
    "dom_ds['frac'].attrs['long_name'] = 'fraction of grid cell that is active'\n",
    "dom_ds['frac'].attrs['units'] = '1'\n",
    "    \n",
    "# Save the output domain to a temporary file.\n",
    "dom_ds.to_netcdf('temp.nc')\n",
    "dom_ds.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cdo gridarea:   0%\b\b\b\b  0%\b\b\b\b  1%\b\b\b\b  2%\b\b\b\b  3%\b\b\b\b  4%\b\b\b\b  5%\b\b\b\b  6%\b\b\b\b  7%\b\b\b\b  8%\b\b\b\b  9%\b\b\b\b 10%\b\b\b\b 11%\b\b\b\b 12%\b\b\b\b 13%\b\b\b\b 14%\b\b\b\b 15%\b\b\b\b 16%\b\b\b\b 17%\b\b\b\b 18%\b\b\b\b 19%\b\b\b\b 20%\b\b\b\b 21%\b\b\b\b 22%\b\b\b\b 23%\b\b\b\b 24%\b\b\b\b 25%\b\b\b\b 26%\b\b\b\b 27%\b\b\b\b 28%\b\b\b\b 29%\b\b\b\b 30%\b\b\b\b 31%\b\b\b\b 32%\b\b\b\b 33%\b\b\b\b 34%\b\b\b\b 35%\b\b\b\b 36%\b\b\b\b 37%\b\b\b\b 38%\b\b\b\b 39%\b\b\b\b 40%\b\b\b\b 41%\b\b\b\b 42%\b\b\b\b 43%\b\b\b\b 44%\b\b\b\b 45%\b\b\b\b 46%\b\b\b\b 47%\b\b\b\b 48%\b\b\b\b 49%\b\b\b\b 50%\b\b\b\b 51%\b\b\b\b 52%\b\b\b\b 53%\b\b\b\b 54%\b\b\b\b 55%\b\b\b\b 56%\b\b\b\b 57%\b\b\b\b 58%\b\b\b\b 59%\b\b\b\b 60%\b\b\b\b 61%\b\b\b\b 62%\b\b\b\b 63%\b\b\b\b 64%\b\b\b\b 65%\b\b\b\b 66%\b\b\b\b 67%\b\b\b\b 68%\b\b\b\b 69%\b\b\b\b 70%\b\b\b\b 71%\b\b\b\b 72%\b\b\b\b 73%\b\b\b\b 74%\b\b\b\b 75%\b\b\b\b 76%\b\b\b\b 77%\b\b\b\b 78%\b\b\b\b 79%\b\b\b\b 80%\b\b\b\b 81%\b\b\b\b 82%\b\b\b\b 83%\b\b\b\b 84%\b\b\b\b 85%\b\b\b\b 86%\b\b\b\b 87%\b\b\b\b 88%\b\b\b\b 89%\b\b\b\b 90%\b\b\b\b 91%\b\b\b\b 92%\b\b\b\b 93%\b\b\b\b 94%\b\b\b\b 95%\b\b\b\b 96%\b\b\b\b 97%\b\b\b\b 98%\b\b\b\b 99%\b\b\b\b100%\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \bcdo gridarea: Processed 3 variables ( 0.63s )\n"
     ]
    }
   ],
   "source": [
    "# This shell command uses cdo step calculates the grid cell area\n",
    "!cdo -O gridarea temp.nc area.nc\n",
    "!rm temp.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:  (lat: 614, lon: 928)\n",
      "Coordinates:\n",
      "  * lon      (lon) float64 -125.0 -124.9 -124.8 -124.8 -124.7 -124.7 -124.6 ...\n",
      "  * lat      (lat) float64 14.66 14.72 14.78 14.84 14.91 14.97 15.03 15.09 ...\n",
      "Data variables:\n",
      "    elev     (lat, lon) float64 nan nan nan nan nan nan nan nan nan nan nan ...\n",
      "    mask     (lat, lon) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...\n",
      "    frac     (lat, lon) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    area     (lat, lon) float64 4.673e+07 4.673e+07 4.673e+07 4.673e+07 ...\n",
      "Attributes:\n",
      "    title: VIC domain data\n",
      "    Conventions: CF-1.6\n",
      "    history: created by jhamman, 2016-06-08 15:01:23.814110\n",
      "    user_comment: VIC domain data\n",
      "    source: generated from VIC North American 1/16 deg. model parameters, see Livneh et al. (2015) for more information\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x114f5c748>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEPCAYAAABIut/fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X28HGV99/HPN/hQLApEiQghoAlFTRXQFi3WggVNfChY\npd6kVbG1rVbxCfWFWu+mtPZurYpPtNUiUrDa2GoVyq0ELYh3i2AUEiEkEGIIOUISS0C0oA05v/uP\nmZMMm9k9s7szuzOz3/frNa+zMzsze82e3d/89pprrksRgZmZtc+ccRfAzMyq4QBvZtZSDvBmZi3l\nAG9m1lIO8GZmLeUAb2bWUg7wZmYlkXSBpG2SvtdjnY9J2iBptaRjqiyPA7yZWXkuBJZ0e1LSC4GF\nEXEk8DrgE1UWxgHezKwkEfEfwD09VjkVuDhd9zpgf0mPr6o8DvBmZqNzKLAlM/+DdFklHODNzEZH\nOcsq6y/mYVXtuAyS3FGOmRUWEXkBtLAjDnt4bJ56sOjq2yLi4D5fYgo4LDM/H7izz30UVusAD7Dr\nriMr2e85H7yb5e94bCX7rlJTyr3kkKNZeecalhxyNAAbYy233fXgXut0s/LONZWWrx91e8+z71v2\nfep8PzfGWhZq8UOW1el97WbQ93ufJ2wY+rU3Tz3IzrsWFlr34U/Y2K3uXORn6gCXAm8EPi/p2cC9\nEbGt74IW5Coaq8RMIFl555rdjzsDUK9g0yv4T7rOoD4zWTl2xXShKY+kzwHXAL8g6Q5JvyvpdZL+\nECAivgJsknQb8EngDVUeiwO8VWIm4HTLNnsts9kN+r71cyKY1JPGNFFoyhMRvx0Rh0TEIyNiQURc\nGBGfjIi/z6xzZkQsioijI+L6Ko9lYgP8CcfvO+4iDKQp5c5m8EsOOZpPffG+nuvmBay6BJi6vuez\nBfkDOajS/Vdl3O/3zthVaGqC2tfBV+XE4x817iIMpO3l7ladM05Nfc/nat64izCQcb/f3bLzJprY\nDN5Gq9+A7aqbYrr9+sl7frZ1LbGLKDQ1wcRm8DY6M9U0M0G+W5Bx8Blct18+dfol1BRtyuAd4G0k\nujXnmy37tOFlm6va7Ha1aJxqV9HYyPVbVeDgVFzn+5qtnrFipgtOTeAM3kZm5uanGaMIOrNVC7XR\nuI+16e95U+rXi1DU+OeIpKjqTlZrjlEHjM4TkTXDPk/YMHRXBZJi45ZivQ8sPGzr0K9XNVfR2FjN\nVv0yjuqZSQ7uw7zfbalK24UKTU1QeYCXdLukNZJukPTtdNmBkq6QdIuklZL2r7ocVl9tCQxtMMzJ\nbZBt6/i/n45iUxOMIoOfBk6MiGMj4rh02buAr0fEUcCVwLtHUA6rsdk6HpvkrLrN6vh/dQbfH+W8\nzqnARenji4CXjqAcVkOzfcHLDAB1zBatfhzg+xPASkmrJP1+uuzxM11kRsRWGLLTDLMCyj5Z+ITR\nTtOhQlMTjKKZ5PERsVXSQcAVkm6hjxFMzvng3bsfn3D8vmPvp8LKV8ef6VZ/37jmfq6+5oHS99uU\n7LyIkTaTlLQc+Anw+yT18tskHQxcFRFPyVk/TtZpDgBWK/30a2+jU1YzyWtvP6LQus8+4va9Xk/S\nUuAjJLUjF0TE+zueXwB8mqTW4m7glRFR2YhOlVbRSHqUpP3Sxz8PvAC4kWRUk9ekq50BXFJlOWwy\n5FWZuCplMpT5Px60Dl7SHOA8YAmwGFgm6ckdq30Q+IeIOBr4M+CvSit4jqqraB4PfCkdW/VhwGcj\n4gpJ3wH+WdLvAXcAv9VrJ3l9l/hmFOs0qgFF/Lmrl7JP4Lti4Lz3OGBDRGwGkLSCpEHJ+sw6TwXe\nChAR35BUaXJbaYCPiE3AMTnLdwAnF9lHZy95Tb8N2qqRPeH75G/DmB68YuNQYEtmfook6GetBl4O\nfFzSy4D9JB0YEfcM+qK9NKYvmrxA7y+xZfnkPxmqrnLrdpF19bX3s+ba+3ttmrdh50XOdwLnSXoN\n8E3gB8CDnRuVpTEBfkYdR/yxenBgb7fZbobb5wnlvE63KpqnPWs/nvas/XbPX/zRHZ2rTAELMvPz\ngYdcQI2Iu0gy+Jnrki+PiB8PXeguGhfgZ2T7uPYX27p9Bib1l16bjrtbYK/q+KYHbya5Clgk6XDg\nLuB0YFl2BUmPBXZE0nzx3SQtairTyADvwG79GMXnZbbsctRGcT1inCeRKl/3f2KwsBgRuySdCVzB\nnmaS6ySdA6yKiMuAE4G/lDRNUkXzxnJKna/23QWfrNMAB3Orh5mglv07m7Z+dvsN8IOcEIqO/gXl\ntYP/8sZi1b8vXbim9t0F1z7Az/QHX7cMySbToNd+/BntXz/BHcoL8F+8ba+Gf7levmh17QN8Y6po\n/AWxuug2vmwvbaoTH4VxVsPuatEwGY0J8Gbj1BlwBsnkHeT7M673anrwG51qp/YB3l8KG7fOYF7m\nqEf+bOcb5/viDN6sxXoF8DLvv3Bw31sdErqdsc9YX79MDvBmGf0E8LxA1M/2owhmdQiYRdWlrEP0\nRVM7jTgS37VqozDK4D7MNv2oQ8Asok7f8WlUaGoCZ/Bm9HcRtTNoDhOcmhKAq1SXzH1GmzJ4B3ib\neEUDTLd1it7wNMxrt1ndjt8XWc1apFtG3k/g6ZX599pPnaomLNGU8VaLcIA36zBMRtnvrftWP87g\nx8A/Za1t/HmupzY1k2zPqcqsRcoaS3acvxKa+gtlOuYUmpqgMRm82SQa9s7Xcf5KaOovlG4jOhUh\naSnwEfZ0F/z+jucPAy4CDkjXeXdEfHXw0vbmAG82QkV7SMzrkrhbJ2dNDaR1NWh2LmkOcB5wEslI\nTqskXRIR2UG33wt8PiI+KekpwFeAJw5Z5K5GEuDTA/8usCUiTpH0D8CvAT8iGbPwNRHxvVGUxWyU\n8jLwblUXect7jWbkgeirMUQ7+OOADRGxGUDSCuBUIBvgp4HHpI8PIBmTtTKjqkh6C7A2Mx/A2yPi\n2Ih4hoO7tcmg9eezBejZ9tnUOu+6GeJO1kOBLZn5qXRZ1jnAqyRtAS4D3lTJQaQqz+AlzQdeBPwF\ncFbmqWZcpTDrQ7a1V7eMfLaqmey6edsXfX0bTLcM/vZVP2Tzqv/qtWle1O8cUWkZcGFEfFjSs4F/\nBBYPUs4iRhFkPwy8k70P9H2SVkv6kKSHj6AcZpXrdbdrr+fLfC1n8sPZGfvkTof+0sEc/0e/uHvK\nMQUsyMzPJ6mLz3ot8M8AEXEt8HOSHlfFcUDFGbykFwPbImK1pBMzT70rIralgf184GzgfXn72Bh7\nana+cc39nHj8oyossVk58i6COrMu1zeuuZ+rr3mg9P0OcSfrKmCRpMOBu4DTSTL2rM3AycBF6UXW\nR0ZEz58Fw6i6iuY5wCmSXgTsCzxa0sUR8WqAiNgp6ULg7d12sFDJr5fky+HgbvWXN0BIFcF9tuqg\ntjvx+Ec9JOH7sw/tKGW/0wNWbETELklnAlewp5nkOknnAKsi4jLgHcD5kt5GcsH1jFIK3UWlAT4i\n3gO8B0DSCSQXVl8t6eCI2CpJwEuBm6osh9k4DRrch+nEbNIk79OGUva1a4i+aCLicuCojmXLM4/X\nAb868Av0aVzt4D+b1jsJWA28fkzlMCtVmQF5kP7pJ+2EUMXxurOxAUTE1cDV6eOTRvW6Zm2V98tg\nUur5qzyRNaUbgiJ8J6tZhQYNuJOWifej6vdmmK4K6sYB3qwkVQ7IPWyfNN106wahLr8ExnGie3C6\nPb1JOsCb1Ui3bgfKvOBa5AaqOgX5UV9sbsp4q0U4wJvVSFXBbNABwUcd5Ae5e7dsw7SiqRsHeLMS\nVFk9M4wq+sSpSl2uO/gi64jV5aeiWVPV/TtUl+AO7Wom2Z5TlVkNVRG4BhkM3IobojfJ2mlEBm9m\nw90RW2d1yt6hXRl87QN83T+cZtYubiY5QnVqrmVWN3m9VlbVZr4qdetzpynVL0XUPsCbWX6ik9d2\nvVPdgzu4iqZKvshqVnODtg1vQnCvo+lQoSmPpKWS1ku6VdLZOc+fK+kGSddLukVSOX0cd+EM3qwE\n3aoZygiysw3jV/brjdrMe1eXqppBM3hJc4DzgJNIRnJaJemSiNg96HZEnJVZ/0zgmOFK21sjArzr\n4a0JOgNUGZ/ZogGvDd+PugT5IapojgM2RMRmAEkrgFOB9V3WXwb8yaAvVkQjArzZJHI1zHg8OPid\nrIcCWzLzUyRBfy+SFgBHAFcO+mJFOMCblaiqQbW7dULWFnXI3GcMkcHnbRhd1j0d+EJEdHu+FA7w\nY+JqJ+tH2z8rM8G9DkG+W4C/Z/UW7lk91WvTKWBBZn4+SV18ntOBNwxSvn44wI9J9sKSWVtlf3nk\nPe5Uh0y+W4Df/+gF7H/0nvh9+8XXdq6yClgk6XDgLpIgvqxzJUlHAQdExF47KJubSY7ZuD/MZlVY\ncsjRe1UrFRk3tg7fhwgVmvbeLnYBZwJXAGuBFRGxTtI5kl6SWfV0YMUojsUZ/Bg5e7e2mu2mrDob\n5k7WiLgcOKpj2fKO+XMGfoE+jSTAp+1DvwNMRcQpko4gOYMdCFwPvCoiHhxFWcxsdJoU2Ge06U7W\nUWXwbwFuBh6Tzr8f+FBE/IukvwNeC3yy28bOdM3qr4nBPM+u6fbUXFd+JJLmAy8CPpVZ/OvAF9PH\nFwG/WXU5zKw6bQnuMHgdfB2NIoP/MPBOYH8ASY8F7omI6fT5KeCQEZTDzGxWrqIpSNKLgW0RsVrS\niTOL2fuGgK6N/TfGWhY94WG86u3bOOH4fTnx+EdVVFoza5IdsZ17+GHp+6321qPRqjqDfw5wiqQX\nAfsCjwY+AuwvaU6axfe6GYCFWgzANecuZvk7XBdvZom5msdc5u2e3xTrStlvm/qDr7QOPiLeExEL\nIuJJJG0/r4yIVwJXAb+VrnYGcEmR/bWpns/M6qlNdfDjulz8LuAsSbcCc4ELxlQOMxtS2xKvYfqD\nr5uR3egUEVcDV6ePNwHPGtVrm1k12hbcAaanmxG8i/CdrGY2kDYGd6Ax1S9FtKdFv5lVrq1BPctV\nNGPkHhjNxqvtQb5NzSQbmcG3/QNmNm7Z3iCzJiG5ciuaGnCQN6vOyjvXTEQwzzNMgJe0VNJ6SbdK\nOrvLOq+QtFbSjZL+scpjaVwVjZlZlQatoUl7zT0POInk5s1Vki6JiPWZdRYBZwO/EhH3SXrc0AXu\nobEZPDiLN6vSpH6/YlqFphzHARsiYnNE7CTpEv3UjnX+APibiLgPICL+q8pjaXSAn9SfkGbjMglB\nf4gqmkOBLZn5qXRZ1i8AR0n6D0nXSFpS0WEADa+iyX7YHOzNytXkUZmGMUQrmryo37m3hwGLgF8j\nGaD7/0laPJPRl63RAd7MqjWJzZK7XUB9YO33eeDmTb02nSIJ2jPyOlKcAr6VdrR4u6RbgCOB7w5c\n4B4c4M1sVpOSvQPQJcDv+9SF7PvUhbvn7/3iVZ2rrAIWSTocuIukg8VlHet8OV12cXqB9Ujg+6WU\nO4cDvJnlmsTsHQavoomIXZLOBK4gub55QUSsk3QOsCoiLouIlZJeIGkt8CDwjoi4p6yyd2pNgJ/U\nD6NZVVbeuWayMvcZQ9zJGhGXA0d1LFveMf924O2Dv0pxrQnwZmZl6NIEspEa3UyyU7fbq81sMJP4\nq7hNXRW0LoOfxA+kmZXInY3VlzN4s3JNXtKkglP9tS6DNzMbijP4enMWb1aOifwuRcGpAVoZ4GFC\nP5hmJZupnpmkapohOhurnUoDvKRHSrpO0g1p38fL0+UXSvp+uvx6SU+vshxmNpyJSphqlsFLembO\nst8osm2lAT4ifgY8LyKOBY4BXijpWenT74iIYyPiGRHxvSpef6I+lGZWjlCxaXTOl/S0mRlJy4D3\nFtmw8iqaiLg/ffhIkou60+n8SN4hB3mz4Uzad0hRbBqh04CLJD1F0h8AbwBeUGTDygO8pDmSbgC2\nAl+LiFXpU++TtFrShyQ9vOpymFl/sjcOTlIdfN2qaCLi+yQdl32RJNi/ICJ+VGRbxYiGEJf0GOBL\nwJuAuyNiWxrYzwdui4j35WwTT+Qpu+cP5CDmat5Arz9RH1CzEtU1g98R27mHH+6e38Q6YshbTCXF\n4Z/4QKF1N7/+nUO/3ixluZGHnkrmAT8CfgYQEbNeuxxZO/h0/MGrgaURcW66bKekC+nR8c5CLR5V\nEc2sQ12DO8BczWMuexK+TbGunB3XpwnkS4bdQdWtaB4naf/08b7AycB6SQenywS8FLipynKYWf/q\nHNwrNV1wqlg6tutmkkR8a/r4iSTjvBaqoikU4CW9v8iyHE8ArpK0GrgOWBkRXwE+K2kNsAZ4LLBX\n9YyZjdZMnfvEd9o3RCsaSUslrZd0q6Szc54/Q9L2tHn49ZJ+r0CJvgjskrQI+HvgMOBzRQ6laBXN\n84HOwr4wZ9lDRMSNwDNylp9U8HVL0+0D67p5swnO1nMM2kJG0hzgPOAkkqH6Vkm6JCLWd6y6IiLe\n3MeupyPiQUkvAz4eER9PG67MqmcGL+mP0or+oyR9LzNtAippuz5q/mDbpPN3oMPgrWiOAzakVSs7\ngRUk1Smd+r0wuzNt+/5q4LJ0WaGWh7Nl8J8Dvgr8JfCuzPIfR8SOPgs5dnnZ+szPUWfyNokc3Et1\nKLAlMz9FEvQ7vUzSc4FbgbMiYmqW/f4u8HrgLyJik6QnAv9YpEA9A3za1vJHpAPHSpoH/Bywn6T9\nIuKOIi9SF70+zNnnHOytzSaybXsfulXRPLDhNn5628aem+Ys69zbpcDn0haErwMuIqnS6Soibgbe\nnJnfBPxVr21mFKqDT/s9OBc4BNgOHA6sA1rZhtEZvbVRZ4Lj7L2LLhdQ9110JPsuOnL3/L2Xf61z\nlSlgQWZ+Pkld/J5dP3SA7fOBWRurSDqSpBblqSQJ9sy+njTbtkWbSb4PeDZwa0Q8keSMc23BbRvJ\nH36zCTV4M8lVwCJJh0t6BMndp5dmV5hpIp46Fbi5QIkuBP4OeBB4HnAxBatoigb4nRFxNzBH0pyI\nuAr4pYLbNsrKO9fszt4nvrmYtYY/x8UN2hdNROwCzgSuANaStJZZJ+kcSTM3Lb1Z0k1pK5gzgdcU\nKNK+EfHvJD0PbI6IPwVeXORYijaTvFfSfsA3Sdqwbwf+u+C2jeIvgtmEG+JO1oi4HDiqY9nyzOP3\nAO/pc7c/TZtgbpB0JvADYL8iGxbN4E8FHgDeBlwObAQK9UfcBg761lT+FTqAmnU2BrwVeBTJhdZn\nAq8kaTI5q0IZfERks/WL+i1dk6y8c42/EGYTbMRdARcRwGdIGrfMtH8/HxiuszFJPyb/XCUgIuIx\n/ZWz/hzcrS38WR7QaAfzKOKzwDuBG+mzF5zZ2sE/eohCtYqbTlpTOLAPqX4Z/A8j4tLZV9vbyLoL\nNrPqObgPTyPoKbJPyyV9Cvh30r7gASLiX2fb0AG+D87irc4c3MtRwzr43wWeTFL/PnP6CcABvmwO\n8mYtV78A/8sRcdTsq+2t8jFZzax6zt5LVL9mktdIeuogGzqDr4g7dLJRcXAvVw2raJ4NrE67af8Z\ne1ox1mdM1kmS/cK5l0ozG9LSQTd0Fc0ABs2YnGlZ2fyZqkDNqmhmxmbtnIps6wx+xFx1Y1ZvNWwm\nOTBn8APqlTnNFrwd3M1qrGYZ/DAqDfCSHinpOkk3SLpR0vJ0+RGSrpV0i6R/kuRfEmZWC4N2Fwwg\naamk9ZJulXR219eQTpM0LekZVR0HVBzgI+JnwPMi4ljgGOCFkp5FMorJh9K2nfcCr62yHFVxXbyB\n/5+tM2AGn3bpex6whGS0u2WSnpyz3n7AmxjBoEmVV9FExP3pw0eS1PkHyagkX0yXXwT8ZtXlMDMr\nYogM/jhgQ3oRdCewgqSr9U5/TpLk/iznuVJVXjWSntW+CywE/oakL/l7I2LmUsYUyVivjTTona2+\nI7b+Ov9HeU1es+vMlsn7/90Qg9evHwpsycxPkQT93SQdA8yPiK9IeufAr1RQ5QE+DeTHSnoM8CXg\nKXmrddt+Y6zd/fhADmKu5pVeRrM8nQE5O1ZAt3sdevFJvVw7Yjv38MPS9ztEK5q8foZ3xzZJAj4M\nnDHLNqUZ2cXNiLhP0tUkd2UdkI7tOk3OyONZC7V4VEUsVdGszppnmEFhugV53xDXv7max1z2JHyb\nYl05O+6Sbv73Hbfx33fc1mvLKWBBZr4ztj2apG7+G2mwPxi4RNIpEXH9MEXupupWNI+TtH/6eF/g\nZJJRxK8Cfitd7QzgkirLUbVBm0w6+FdrZri6st7nvOy9zP12m+/FJ4MKdLmo+vOHLWLec5bunnKs\nAhZJOlzSI4DTgd39uEfEfRExLyKeFBFPJLnI+htVBXeoPoN/AnBRWg8/B/h8Wve0Dlgh6c+BG4AL\nKi5H5Xr9/PYwgNUqUg/emSHPbFO02qTM/58/C/U2aF80EbErHRT7CpJ4d0FErJN0DrAqIi7r3ISK\nq2gUUd8W+5LiZJ027mL0pVuwyAaUottYf4YJnEX/B6MMznUsU519Pb5AxHDj7UmKxe84t9C6az94\n1tCvVzXfyVqyblUC/WSL1r+yM+xsdcy46seLHpM/U+Ua5kanunGAH6G8DN5fzuGUVcee93/Jaykz\n6mzZ2fkYuKsCm02vTN7KUcXFzl6BfFzBtsjr+nNVnjZl8O4DZszcPG4wzmytMg0J3kU4g69YP4HI\nQau7bnXitocThHK0KYN3gB+DlXeu6dnaxvbmpqY2Mi2qg3cVTYkchKox6e+pM/MRa0jwLsIZfA1N\nekAzGydX0ViuMofjK/s2+6aa9OPvl7P9EriKxmZT5IvmKp3e/N4kfIPcaHlMVuupn8zbX1yzenEV\njZXKQX5vzt5tbFpUReMAX6Gyu6m1yeUuhEfIAd6KKlpdM9uX0hdcrRd/NsozTBWNpKWS1ku6VdLZ\nOc+/TtL3JN0g6Zt5g3KXyQG+Rpx5WS8O4iMyYAafjntxHrCEZOSmZTkB/LMR8fSIOBb4AMkQfpVx\ngB+RIs0ei7SWmJRM3ie7/vTzfs3cSZ3dxu/3HoooNOU4DtgQEZsjYiewAjg1u0JE/CQzux9QaZsd\nN5Mcg27t5fv5krnpnA2qW0+ZbrabGKKZ5KHAlsz8FEnQf+j+pTcAZwEPB3594FcrwBn8GPnLZP2a\nlF9wYzX4Rda80Z32WjMi/jYiFgFnA/+7lDJ34Qx+zPIycWdSNmr+NbhHtwuo9229jR9v3dhr0ylg\nQWZ+PnBnj/U/D3yiz+L1xRl8DTiY27j5M5jRJWN/zOMXcejRS3ZPOVYBiyQdLukRwOnApdkVJC3K\nzL4EuLWKQ5hRaYCXNF/SlZJulnSjpDely5dLmpJ0fTotrbIcZm1TRUB2kE8M2kwyInYBZwJXAGuB\nFRGxTtI5kl6SrnampJskXQ+8FTijymOpuormQeCsiFgtaT/gu5K+lj53bkQUG758AhWppvGF1snW\n7f/vQD2kIW5iiojLgaM6li3PPH7r4HvvX6UZfERsjYjV6eOfAOtIrjRD/gUJ61Nbv8xtPa6y+X0q\nn/uiGYCkI4BjgOvSRW+UtFrSpyTtP6pyNEnR7NxfcrPyaDoKTU0wkgCfVs98AXhLmsn/LbAwIo4B\ntgJdq2o2xtrd047YPoriNlKbms+15ThGaZD3rOnVezti+0PiQ2la1BeNIv+OrPJeQHoYcBnw1Yj4\naM7zhwP/FhFPz3kuTtZplZavLnp92fr98jb5i+vgboP6enyBiBiq6ldS/MorPlho3W/98zuGfr2q\njSKD/zRwcza4Szo48/zLgJtGUI5ac2Azq4kWZfCVtqKR9Bzgd4AbJd1A8ra8B/htSceQ9MNwO/C6\nKsvRFHmtIgYJ/GUOHWg2aZpyAbWISgN8RPwnsE/OU5dX+bpN1hnkh7mr1c0ozQZQcbX1KPlO1hrK\nu2DqQG02GpouNjWBA3wDzAT3QYK86/bN+uN28DYSZWXyDvJmfYgoNjWAA/yEcJA3K8YZvI2U6+PN\nRqhFzSQd4BuijAzcWbzZ7JzBW2M5yJvNYjqKTQ3gAN9AwwZpB3mz7oZpJilpqaT1km6VdHbO82+T\ntDbtaPFrkg6r8lg8ZF+DODCbjcCALWQkzQHOA04iGapvlaRLImJ9ZrXrgWdGxE8lvR74AMnIT5Vw\nBj+hfLIwyzdEHfxxwIaI2BwRO4EVwKnZFSLi6oj4aTp7LXvGx6iEA/wEc5A3yzF4K5pDgS2Z+Sl6\nB/DXAl8durw9uIrGzCxDXapo7tnxfe695/s9N81ZlrszSa8Engmc0G/5+uEAb2aW1eUC6oEHPIkD\nD3jS7vnbN/175ypTwILM/HySuviHkHQy8G7g19KqnMo4wJuZZXTL4AtYBSxKBzG6i+Ti6bKH7Fs6\nFvgEsCQi7h6mnEW4Dn7C1a0e3nfp2tgN2A4+InYBZwJXAGuBFRGxTtI5kl6SrvbXwM8D/yLpBklf\nrvJQnMGbmWUMc5dqRFwOHNWxbHnm8fMH33v/HODNA4MMoNv7NfNe1u2XkfWhIT1FFuEAb7VTxwBZ\n9ASY7bu/bsdgxTRlMI8iHODNevAvmwnkDN7axEEs3zDvi7P3BmtPfK+2FY2k+ZKulHSzpBslvTld\nfqCkKyTdImmlpP2rLIc1z7hPOuN+fRsfRRSamqDqZpIPAmdFxFOBXwHeKOnJwLuAr0fEUcCVJI3+\nbUycbZbL72fD7YpiUwNUGuAjYmtErE4f/wRYR3J316nARelqFwEvrbIc1t3KO9fUNluta7l6mQnu\nTSy7JdqUwY+sDl7SEcAxJD2oPT4itkFyEpB00KjKYTYKzuIbrCHBu4iR3MkqaT/gC8Bb0ky+Pe9g\nw9U9EDkTtpGLKDY1QOUZvKSHkQT3z0TEJenibZIeHxHbJB0MbO+2/cZYu/vxgRzEXM2rtLxWP25T\nbnl2xHbu4Yfl77hF7eBHkcF/Grg5Ij6aWXYp8Jr08RnAJZ0bzVioxbsnB/dqNCF4jjKTH/S1mvA+\ntslczXtIfChLm+rgq24m+Rzgd4BfTzvWuV7SUuD9wPMl3QKcDPxVleWw2Tk4DcfvX4u4iqaYiPhP\nYJ8uT5+t+9KYAAAJI0lEQVRc5WubjZKrkVpkuj11NO4u2Hare4AaRTXNMK/hC8ItMV1wyiFpqaT1\nkm6VdHbO88+V9F1JOyW9rKIj2M0B3h5iySFH1z7QV2XYAD2p71vbDFoHL2kOcB6wBFgMLEtv7Mza\nTHLd8bNVHwc4wFsXkxasnH3bboPXwR8HbIiIzelQfCtIburM7DruiIibGFFTcXc2Zl3V8a7MzrKM\n+0Tk/t9bKGe0poIOBbZk5qdIgv7YOMDbRMmeIByULVeXFjJ3338HOx64o9eWyttbGUUalAO8zarO\nIz7NlKtbsO5V7rLq3H2iaJkuAf6x+x7GY/c9bPf8xh3XdK4yBSzIzM8H7iy5dH1xgLdC6hzk4aHN\nFEdRTgf1Fts1cDPJVcAiSYcDdwGnA8t6rJ+X8ZfKAd4Kywa1Ogb7OpbJGigGC/ARsUvSmcAVJA1Y\nLoiIdZLOAVZFxGWSfgn4EnAA8BJJfxoRTyur6J0c4G0gnRlsG4Krs3IDhrpLNSIuB47qWLY88/g7\nwGGd21XFzSStFA6O1hrTUWxqAAd4K42DvLVCi/qicYC3UjnIW+O1KMC7Dt5KV8cbpGbjE5PttmvX\nuEtQGmfwVpkm9WvTpJORVaxFGbwDvFWuKYHeQd4AB3izQTQh0DvIW5ta0bgO3kauCTdM1f1EZNWJ\nAW90qiMHeBurMi7INvGirtVYQ7LzIlxFY7UwaMac3c5Zt5WiRXXwzuCtNvrp0KxbMO+3+scnBdtL\ni5pJVhrgJV0AvATYFhFPT5ctB/4A2J6u9p60/wazWfu46Scgz3bCcHC3PNGiQberzuAvBD4OXNyx\n/NyIOLfi17YWGDYI172bY6uhhlS/FFFpHXxE/AdwT85TlfeDbDYj7yTh7N26GqKZpKSlktZLulXS\n2TnPP0LSCkkbJH1L0oK8/ZRlXBdZ3yhptaRPSdp/HAXYEdtnX6mGXO7BDBPQx132QbncA4rpYlMH\nSXOA84AlwGJgmaQnd6z2WmBHRBwJfAT46yoPZRwB/m+BhRFxDLAVGEtVzT38cBwvOzSXe3CDtr+v\nQ9kH4XIPJqaj0JTjOGBDRGyOiJ3ACuDUjnVOBS5KH38BOKmyA2EMrWgiIvvfOx/4t17rb4y1ux8f\nyEHM1byKSmaTxFU0zbcjtldzMhj8RqdDgS2Z+SmSoJ+7TjoC1L2S5kbEjkFftJdRBHiRqXOXdHBE\nbE1nXwbc1GvjhVpcYdFs0viia3vM1Tzmsifh2xTrStlvDN5MMu/aYmeq37mOctYpjaLCK8aSPgec\nCDwW2AYsB54HHANMA7cDr4uIbV22b8/lbDOrXEQM1YBD0u3A4QVX3xYRB2e2fTbwpxGxNJ1/V1Kk\neH9mna+m61wnaR/groiorFqi0gBvZjYp0oB9C0m9+l3At4FlEXt+Wkh6A/CLEfEGSacDL42I06sq\nk+9kNTMrQVqnfiZwBUkDlgsiYp2kc4BVEXEZcAHwGUkbgLuByoI7OIM3M2ut1nc2Juk0STdJ2iXp\nGZnlJ0v6jqQ1klZJel7muWdI+l56s8JHxlPynmWfK+lKST+W9LGObcZe9m7lTp97d3qTxzpJL8gs\n73mDyKhJerqka9LPxyWS9ss8l3sMdSDp6PQGmhskfVvSL2ee+1ha7tWSjhlnOfOkNwBdn06bJF2f\nea6273mtRUSrJ+Ao4EjgSuAZmeVHAwenjxcDU5nnrgOOSx9/BVhSs7I/Cjge+EPgYx3bjL3sPcr9\nFOAGkqrBI4DbSFoRzEkfHw48HFgNPHnMn5tvA7+aPn4N8Gfp46fmHcM4y9pR7pXAC9LHLwSuSh+/\nCPi/6eNnAdeOu6yzHMcHgff2+tyMu4xNmFqfwUfELRGxgY7mSRGxJtLmmhGxFnikpIdLOhh4dER8\nO131YuClIy30njJ2K/v9EXEN8LPs8rqUvVu5SW7yWBERD0bE7cAGknbCRW4QGbVfiKSrDYCvAy9P\nH59C/jHUxTQwc3f4AcAP0senkPYJFRHXAftLevzoi1fYK4DPpY+7fW5sFq0P8EVIOg24IQ0uh5Lc\noDBjKl3WBHUve+eNID9Il+XdIDLuct8k6TfSx68A5qePux1DXbwN+KCkO0hug393urzu5d5N0nOB\nrRHx/XRRY8peN61oRSPpa0A2G5m5eeCPI6LnnbKSFgN/CTw/s22n6m5EGKLsebvLWVZJ2Qcsd7fy\n5SUalV/973UMwO8BH5f0J8ClwP9k1uk00pYKs5T7ZOAtEfHlNHH5NMlne+zlhsKfm2XAP3Ws08mt\nQwpoRYCPiOfPvtbeJM0H/hV4VfrTD5Ls8bDMavOBO4cqYA+Dlr2LkZV9wHJ3K5+ABTnLK1XgGJYA\nSDoSeHG6bKSfjzy9yi3pMxHxlnS9L0j6VPrU2Mudlqnne562JX8ZkL04X4uyN9GkVdFku0zYH7gM\neFdEXDuzPK2Xv0/ScZIEvBq4ZOQl3Vu3O/R2L69p2bPlvhQ4Pe0y9YnAIpKLmauARZIOl/QIkrbB\nl46+qHtIOij9Owd4L/CJ9Klux1AXP5B0AoCkk0jqqyEp96vT5c8G7o0ud5CP2fOBdRGRDeB1f8/r\na9xXeaueSC4ybgEeILm77Kvp8j8GfgxcT3KF/nrgcelzzwRuJPlyfLRuZU+f2wT8F3AfcAdpq5M6\nlH2Wcr+bpBXEOtLWHunypSR3AW4gOemO+3Pz5rQ864H/0/Fc7jHUYSJpXfWd9DP9LeDYzHPnpeVe\nQ6Z1U50mkkGC/jBneW3f8zpPvtHJzKylJq2KxsxsYjjAm5m1lAO8mVlLOcCbmbWUA7yZWUs5wJuZ\ntZQDvNWSpB+PuwxmTecAb3XlGzTMhuQAb7Un6QOSbkwH33hFuuwESVdJ+pd0EIjPjLucZnXTis7G\nrL0kvRx4ekQ8TdI8YJWkq9OnjyEZgGMr8J+Sjo+kn3wzwxm81d9zSLuOjYjtwDeAmWHovh0Rd0XS\n38ZqktF+zCzlAG9119mLZnY+O6LVLvyL1OwhHOCtrmYC+TeB/yVpTtqF73NxV7FmhTjjsboKgIj4\nUtp/+RqS8UbfGRHbJT0lb30z28PdBZuZtZSraMzMWsoB3syspRzgzcxaygHezKylHODNzFrKAd7M\nrKUc4M3MWsoB3syspf4/8ELa3v00WXcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1041d5160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This step extracts the area from the temporary area.nc file\n",
    "area = xr.open_dataset('area.nc')['cell_area']\n",
    "dom_ds['area'] = area\n",
    "\n",
    "# Write the final domain file\n",
    "dom_ds.to_netcdf(dom_file)\n",
    "dom_ds.close()\n",
    "\n",
    "# Document the domain and plot\n",
    "print(dom_ds)\n",
    "dom_ds['mask'].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Livneh Parameters\n",
    "\n",
    "VIC 5 uses the same parameters as VIC 4. The following steps will read/parse the ASCII formatted parameter files and construct the netCDF formatted parameter file.  We'll use the domain file constructed in the previous step to help define the spatial grid.\n",
    "\n",
    "Steps:\n",
    "\n",
    "1. Read the soil/snow/veg/veglib files\n",
    "2. Read the target grid (domain file)\n",
    "3. Map the parameters to the spatial grid defined by the domain file\n",
    "4. Write the parameters to a netCDF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading ./Livneh_0.0625_NLDAS/vic.nldas.mexico.soil.txt\n",
      "reading ./Livneh_0.0625_NLDAS/vic.nldas.mexico.snow.txt.L13\n",
      "reading ./Livneh_0.0625_NLDAS/vic.nldas.mexico.veg.txt\n",
      "reading ./Livneh_0.0625_NLDAS/LDAS_veg_lib\n",
      "Reading input data variables:  odict_keys(['lon', 'lat', 'elev', 'mask', 'frac', 'area']) from file: ./domain.vic.global0.0625deg.20160608.nc\n",
      "gridding params now...\n",
      "creating 2d coordinate arrays\n",
      "236213 masked values\n",
      "dont have units or description for elev\n",
      "writing var: run_cell\n",
      "writing var: gridcell\n",
      "writing var: lats\n",
      "writing var: lons\n",
      "writing var: infilt\n",
      "writing var: Ds\n",
      "writing var: Dsmax\n",
      "writing var: Ws\n",
      "writing var: c\n",
      "writing var: expt\n",
      "writing var: Ksat\n",
      "writing var: phi_s\n",
      "writing var: init_moist\n",
      "writing var: elev\n",
      "writing var: depth\n",
      "writing var: avg_T\n",
      "writing var: dp\n",
      "writing var: bubble\n",
      "writing var: quartz\n",
      "writing var: bulk_density\n",
      "writing var: soil_density\n",
      "writing var: off_gmt\n",
      "writing var: Wcr_FRACT\n",
      "writing var: Wpwp_FRACT\n",
      "writing var: rough\n",
      "writing var: snow_rough\n",
      "writing var: annual_prec\n",
      "writing var: resid_moist\n",
      "writing var: fs_active\n",
      "writing var: cellnum\n",
      "writing var: AreaFract\n",
      "writing var: elevation\n",
      "writing var: Pfactor\n",
      "writing var: Nveg (614, 928)\n",
      "writing var: Cv (12, 614, 928)\n",
      "writing var: root_depth (12, 3, 614, 928)\n",
      "writing var: root_fract (12, 3, 614, 928)\n",
      "writing var: LAI (12, 12, 614, 928)\n",
      "writing var: overstory (12, 614, 928)\n",
      "writing var: rarc (12, 614, 928)\n",
      "writing var: rmin (12, 614, 928)\n",
      "writing var: wind_h (12, 614, 928)\n",
      "writing var: RGL (12, 614, 928)\n",
      "writing var: rad_atten (12, 614, 928)\n",
      "writing var: wind_atten (12, 614, 928)\n",
      "writing var: trunk_ratio (12, 614, 928)\n",
      "writing var: albedo (12, 12, 614, 928)\n",
      "writing var: veg_rough (12, 12, 614, 928)\n",
      "writing var: displacement (12, 12, 614, 928)\n"
     ]
    }
   ],
   "source": [
    "soil_file = os.path.join(ldpath, 'vic.nldas.mexico.soil.txt')\n",
    "snow_file = os.path.join(ldpath, 'vic.nldas.mexico.snow.txt.L13')\n",
    "veg_file = os.path.join(ldpath, 'vic.nldas.mexico.veg.txt')\n",
    "vegl_file = os.path.join(ldpath, 'LDAS_veg_lib')\n",
    "out_file = os.path.join(opath, 'livneh_nldas.mexico_vic_5.0.0_parameters.nc')\n",
    "\n",
    "# Set options that define the shape/type of parameters\n",
    "cols = gp.Cols(nlayers=3,\n",
    "               snow_bands=5,\n",
    "               organic_fract=False,\n",
    "               spatial_frost=False,\n",
    "               spatial_snow=False,\n",
    "               july_tavg_supplied=False,\n",
    "               veglib_fcan=False,\n",
    "               veglib_photo=False)\n",
    "n_veg_classes = 11\n",
    "vegparam_lai = True\n",
    "lai_src = 'FROM_VEGPARAM'\n",
    "\n",
    "# ----------------------------------------------------------------- #\n",
    "\n",
    "# Read the soil parameters\n",
    "soil_dict = gp.soil(soil_file, c=gp.Cols(nlayers=3))\n",
    "\n",
    "# Read the snow parameters\n",
    "snow_dict = gp.snow(snow_file, soil_dict, c=cols)\n",
    "\n",
    "# Read the veg parameter file\n",
    "veg_dict = gp.veg(veg_file, soil_dict,\n",
    "                  vegparam_lai=vegparam_lai, lai_src=lai_src, \n",
    "                  veg_classes=n_veg_classes)\n",
    "\n",
    "# Read the veg library file\n",
    "veg_lib, lib_bare_idx = gp.veg_class(vegl_file, c=cols)\n",
    "\n",
    "# Determine the grid shape\n",
    "target_grid, target_attrs = gp.read_netcdf(dom_file)\n",
    "for old, new in [('lon', 'xc'), ('lat', 'yc')]:\n",
    "    target_grid[new] = target_grid.pop(old)\n",
    "    target_attrs[new] = target_attrs.pop(old)\n",
    "\n",
    "# Grid all the parameters\n",
    "grid_dict = gp.grid_params(soil_dict,\n",
    "                           target_grid,\n",
    "                           version_in='4.1.2.c', \n",
    "                           vegparam_lai=vegparam_lai,\n",
    "                           lib_bare_idx=lib_bare_idx, \n",
    "                           lai_src=lai_src,\n",
    "                           veg_dict=veg_dict,\n",
    "                           veglib_dict=veg_lib, \n",
    "                           snow_dict=snow_dict,\n",
    "                           lake_dict=None)\n",
    "\n",
    "# Write a netCDF file with all the parameters\n",
    "gp.write_netcdf(out_file, \n",
    "                target_attrs,\n",
    "                target_grid=target_grid,\n",
    "                vegparam_lai=vegparam_lai,\n",
    "                lai_src=lai_src, \n",
    "                soil_grid=grid_dict['soil_dict'],\n",
    "                snow_grid=grid_dict['snow_dict'],\n",
    "                veg_grid=grid_dict['veg_dict'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global 1/2 deg. Parameters\n",
    "\n",
    "For the global 1/2 deg. parameters, we will follow the same steps as for the Livneh case with one exception. We don't have a domain file this time so we'll use `tonic`'s `calc_grid` function to make one for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading ./Nijssen_0.5_Global/global_soil_param_new\n",
      "reading ./Nijssen_0.5_Global/global_snowbands_new\n",
      "reading ./Nijssen_0.5_Global/global_veg_param_new\n",
      "reading ./Nijssen_0.5_Global/world_veg_lib.txt\n",
      "Calculating grid size now...\n",
      "found 717 unique lons\n",
      "found 278 unique lats\n",
      "assuming grid is meant to be global...\n",
      "creating 2d coordinate arrays\n",
      "Created a target grid based on the lats and lon in the soil parameter file\n",
      "Grid Size: (360, 720)\n",
      "gridding params now...\n",
      "creating 2d coordinate arrays\n",
      "197855 masked values\n",
      "writing var: run_cell\n",
      "writing var: gridcell\n",
      "writing var: lats\n",
      "writing var: lons\n",
      "writing var: infilt\n",
      "writing var: Ds\n",
      "writing var: Dsmax\n",
      "writing var: Ws\n",
      "writing var: c\n",
      "writing var: expt\n",
      "writing var: Ksat\n",
      "writing var: phi_s\n",
      "writing var: init_moist\n",
      "writing var: elev\n",
      "writing var: depth\n",
      "writing var: avg_T\n",
      "writing var: dp\n",
      "writing var: bubble\n",
      "writing var: quartz\n",
      "writing var: bulk_density\n",
      "writing var: soil_density\n",
      "writing var: off_gmt\n",
      "writing var: Wcr_FRACT\n",
      "writing var: Wpwp_FRACT\n",
      "writing var: rough\n",
      "writing var: snow_rough\n",
      "writing var: annual_prec\n",
      "writing var: resid_moist\n",
      "writing var: fs_active\n",
      "writing var: cellnum\n",
      "writing var: AreaFract\n",
      "writing var: elevation\n",
      "writing var: Pfactor\n",
      "writing var: Nveg (360, 720)\n",
      "writing var: Cv (12, 360, 720)\n",
      "writing var: root_depth (12, 2, 360, 720)\n",
      "writing var: root_fract (12, 2, 360, 720)\n",
      "writing var: LAI (12, 12, 360, 720)\n",
      "writing var: overstory (12, 360, 720)\n",
      "writing var: rarc (12, 360, 720)\n",
      "writing var: rmin (12, 360, 720)\n",
      "writing var: wind_h (12, 360, 720)\n",
      "writing var: RGL (12, 360, 720)\n",
      "writing var: rad_atten (12, 360, 720)\n",
      "writing var: wind_atten (12, 360, 720)\n",
      "writing var: trunk_ratio (12, 360, 720)\n",
      "writing var: albedo (12, 12, 360, 720)\n",
      "writing var: veg_rough (12, 12, 360, 720)\n",
      "writing var: displacement (12, 12, 360, 720)\n"
     ]
    }
   ],
   "source": [
    "soil_file = os.path.join(gdpath, 'global_soil_param_new')\n",
    "snow_file = os.path.join(gdpath, 'global_snowbands_new')\n",
    "veg_file = os.path.join(gdpath, 'global_veg_param_new')\n",
    "vegl_file = os.path.join(gdpath, 'world_veg_lib.txt')\n",
    "out_file = os.path.join(gdpath, 'global_0.5deg.vic_5.0.0_parameters.nc')\n",
    "\n",
    "# Set options that define the shape/type of parameters\n",
    "cols = gp.Cols(nlayers=3,\n",
    "               snow_bands=5,\n",
    "               organic_fract=False,\n",
    "               spatial_frost=False,\n",
    "               spatial_snow=False,\n",
    "               july_tavg_supplied=False,\n",
    "               veglib_fcan=False,\n",
    "               veglib_photo=False)\n",
    "n_veg_classes = 11\n",
    "root_zones = 2\n",
    "vegparam_lai = True\n",
    "lai_src = 'FROM_VEGPARAM'\n",
    "\n",
    "# ----------------------------------------------------------------- #\n",
    "\n",
    "# Read the soil parameters\n",
    "soil_dict = gp.soil(soil_file, c=cols)\n",
    "\n",
    "# Read the snow parameters\n",
    "snow_dict = gp.snow(snow_file, soil_dict, c=cols)\n",
    "\n",
    "# Read the veg parameter file\n",
    "veg_dict = gp.veg(veg_file, soil_dict,\n",
    "                  vegparam_lai=vegparam_lai, lai_src=lai_src, \n",
    "                  veg_classes=n_veg_classes, max_roots=root_zones)\n",
    "\n",
    "# Read the veg library file\n",
    "veg_lib, lib_bare_idx = gp.veg_class(vegl_file, c=cols)\n",
    "\n",
    "# Determine the grid shape\n",
    "target_grid, target_attrs = gp.calc_grid(soil_dict['lats'], soil_dict['lons'])\n",
    "\n",
    "# Grid all the parameters\n",
    "grid_dict = gp.grid_params(soil_dict, target_grid, version_in='4', \n",
    "                           vegparam_lai=vegparam_lai, lai_src=lai_src,\n",
    "                           lib_bare_idx=lib_bare_idx, \n",
    "                           veg_dict=veg_dict, veglib_dict=veg_lib, \n",
    "                           snow_dict=snow_dict, lake_dict=None)\n",
    "\n",
    "# Write a netCDF file with all the parameters\n",
    "gp.write_netcdf(out_file,\n",
    "                target_attrs,\n",
    "                target_grid=target_grid,\n",
    "                vegparam_lai=vegparam_lai,\n",
    "                lai_src=lai_src,\n",
    "                soil_grid=grid_dict['soil_dict'],\n",
    "                snow_grid=grid_dict['snow_dict'],\n",
    "                veg_grid=grid_dict['veg_dict'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global Domain File\n",
    "\n",
    "Since the global soil parameters didn't come with a domain file that we could use, we'll construct one using the output from `tonic`'s `calc_grid` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dom_ds = xr.Dataset()\n",
    "\n",
    "# Set global attributes\n",
    "dom_ds.attrs['title'] = 'VIC domain data'\n",
    "dom_ds.attrs['Conventions'] = 'CF-1.6'\n",
    "dom_ds.attrs['history'] = 'created by %s, %s' % (user, now)\n",
    "dom_ds.attrs['user_comment'] = 'VIC domain data'\n",
    "dom_ds.attrs['source'] = 'generated from VIC Global 0.5 deg. model parameters, see Nijssen et al. (2001) for more information'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dom_file = os.path.join(opath, 'domain.vic.global0.5deg.%s.nc' % now.strftime('%Y%m%d'))\n",
    "\n",
    "# Get the mask variable\n",
    "dom_ds['mask'] = xr.DataArray(target_grid['mask'], coords={'lat': target_grid['yc'],\n",
    "                                                                 'lon': target_grid['xc']},\n",
    "                                    dims=('lat', 'lon', ))\n",
    "\n",
    "# For now, the frac variable is going to be just like the mask\n",
    "dom_ds['frac'] = dom_ds['mask'].astype(np.float)\n",
    "dom_ds['frac'].attrs['long_name'] = 'fraction of grid cell that is active'\n",
    "dom_ds['frac'].attrs['units'] = '1'\n",
    "\n",
    "# Set variable attributes\n",
    "for k, v in target_attrs.items():\n",
    "    if k == 'xc':\n",
    "        k = 'lon'\n",
    "    elif k == 'yc':\n",
    "        k = 'lat'\n",
    "    dom_ds[k].attrs = v\n",
    "    \n",
    "# Write temporary file for gridarea calculation\n",
    "dom_ds.to_netcdf('temp.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cdo gridarea:   0%\b\b\b\b  0%\b\b\b\b  1%\b\b\b\b  2%\b\b\b\b  3%\b\b\b\b  4%\b\b\b\b  5%\b\b\b\b  6%\b\b\b\b  7%\b\b\b\b  8%\b\b\b\b  9%\b\b\b\b 10%\b\b\b\b 11%\b\b\b\b 12%\b\b\b\b 13%\b\b\b\b 14%\b\b\b\b 15%\b\b\b\b 16%\b\b\b\b 17%\b\b\b\b 18%\b\b\b\b 19%\b\b\b\b 20%\b\b\b\b 21%\b\b\b\b 22%\b\b\b\b 23%\b\b\b\b 24%\b\b\b\b 25%\b\b\b\b 26%\b\b\b\b 27%\b\b\b\b 28%\b\b\b\b 29%\b\b\b\b 30%\b\b\b\b 31%\b\b\b\b 32%\b\b\b\b 33%\b\b\b\b 34%\b\b\b\b 35%\b\b\b\b 36%\b\b\b\b 37%\b\b\b\b 38%\b\b\b\b 39%\b\b\b\b 40%\b\b\b\b 41%\b\b\b\b 42%\b\b\b\b 43%\b\b\b\b 44%\b\b\b\b 45%\b\b\b\b 46%\b\b\b\b 47%\b\b\b\b 48%\b\b\b\b 49%\b\b\b\b 50%\b\b\b\b 51%\b\b\b\b 52%\b\b\b\b 53%\b\b\b\b 54%\b\b\b\b 55%\b\b\b\b 56%\b\b\b\b 57%\b\b\b\b 58%\b\b\b\b 59%\b\b\b\b 60%\b\b\b\b 61%\b\b\b\b 62%\b\b\b\b 63%\b\b\b\b 64%\b\b\b\b 65%\b\b\b\b 66%\b\b\b\b 67%\b\b\b\b 68%\b\b\b\b 69%\b\b\b\b 70%\b\b\b\b 71%\b\b\b\b 72%\b\b\b\b 73%\b\b\b\b 74%\b\b\b\b 75%\b\b\b\b 76%\b\b\b\b 77%\b\b\b\b 78%\b\b\b\b 79%\b\b\b\b 80%\b\b\b\b 81%\b\b\b\b 82%\b\b\b\b 83%\b\b\b\b 84%\b\b\b\b 85%\b\b\b\b 86%\b\b\b\b 87%\b\b\b\b 88%\b\b\b\b 89%\b\b\b\b 90%\b\b\b\b 91%\b\b\b\b 92%\b\b\b\b 93%\b\b\b\b 94%\b\b\b\b 95%\b\b\b\b 96%\b\b\b\b 97%\b\b\b\b 98%\b\b\b\b 99%\b\b\b\b100%\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \bcdo gridarea: Processed 2 variables ( 0.28s )\n"
     ]
    }
   ],
   "source": [
    "# This step calculates the grid cell area\n",
    "!cdo -O gridarea temp.nc area.nc\n",
    "!rm temp.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:  (lat: 360, lon: 720)\n",
      "Coordinates:\n",
      "  * lat      (lat) float64 -89.75 -89.25 -88.75 -88.25 -87.75 -87.25 -86.75 ...\n",
      "  * lon      (lon) float64 -179.8 -179.2 -178.8 -178.2 -177.8 -177.2 -176.8 ...\n",
      "Data variables:\n",
      "    mask     (lat, lon) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...\n",
      "    frac     (lat, lon) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    area     (lat, lon) float64 nan nan nan nan nan nan nan nan nan nan nan ...\n",
      "Attributes:\n",
      "    title: VIC domain data\n",
      "    Conventions: CF-1.6\n",
      "    history: created by jhamman, 2016-06-08 15:01:23.814110\n",
      "    user_comment: VIC domain data\n",
      "    source: generated from VIC Global 0.5 deg. model parameters, see Nijssen et al. (2001) for more information\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x11d88d7b8>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEPCAYAAACk43iMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X/0XHdd5/HnK5WuQPnRYH9AmgRoSlp6oLVglvUnq9AE\ntxoOVDbxeCg/dN2FrgdRtqDuCXXxYPWIqF2WFaqnCBgUFCMLJfxQWcVCoE1bQtqmaZvmS39Kw48F\n5KTf73v/mJn0dnLvzL1z7525P16Pc+Z8Z+73zszn3pn5vO/ntyICMzPrp1WLToCZmS2Og4CZWY85\nCJiZ9ZiDgJlZjzkImJn1mIOAmVmPOQiYmc2JpKsk3Sfpxgn7/KGkA5L2Sjq/7jQ5CJiZzc+fApuz\n/inpRcCZEXEW8IvAO+tOkIOAmdmcRMQ/Akcm7LIVeM9w388BT5B0Wp1pchAwM2uONcDhxOOvDLfV\nxkHAzKw5lLKt1rl9vqfOF6+bJE98ZGa5RURaJpvbU9c+Kg4tPZR39/si4vSCb7EErE08PgO4u+Br\nFNLqIADwAl1cy+sejH2cqXNree1ZOU35OE35NTFddaXpk/HB0q9xaOkhjt5zZq59H/Xkg1l1+SL9\nih9gF/Ba4AOSngd8LSLuK5zQAlofBMzM5mk5VmZ+rqT3A88HniTpLmAHcCIQEfHHEfFRST8p6Tbg\nW8ArK0jyRA4CZmYFrJSooo+In82xz6Uzv8EMHAQynMwpi07CcZymfJym/JqYriamKeloLC86CZVy\nEMiwWqcuOgnHcZrycZrya2K6mpimpDIlgSZyEDAzK2DZQcDMrL9cEjCb4ON338Dmp5z3iMfTJPc3\na7rljq3L7iDQQaOMd/NTzsuVCWfJypxHGX3a+8zynuOBw6zJZu8g2kwOAh2QlenOGgBGGXLa89P+\nl3U/67WnpXc8IGRt7wIHwPbpWpuAosVFG0lR14jhpktmpFmZ9ijDzXN1Pinjn5dJ6WxrRpkWwJKf\nSdoxV32sk4J5n3wyPlh62ghJcfBwvpkgzlx7b+n3mwcHgRaalFEXqY5JywgWGQTG1ZVRFTl/Xcgs\ns4JOF46tiKqCwK2Hn5xr32esvacVQcDVQQ016Yc6qeF12pV/W374daYzrQ0j+X5tOUd55fm+WH4r\n7b1uTuUg0FBpP9ysq7jxq7w8V/NNuuIfN69MKXk+x6+WuyirraXLx1yH5cy539rJ1UELlKd4nvYD\nndaAmqXp1T9wfECbR+aUfJ++VpX0QVXVQTfelW+Nl2ev+0orqoO8qMwcJa+8kjY/5bzMKpzx7Wn7\n5L36T9O0zC6ZIc+zRJC834RzUuYzneV90t6vaRcITbESynVrC5cE5mxaj5Cs3iTj+43/L28pYtY0\nz1sTMuJ5mEfPnbKfX1XpmVatWbeqSgJfPLR2+o7Ac9YfbkVJwG0Cc5Z1xV+kx8q0QDH+um3LULte\nNz9SZy+kJl/Fj3+X29Zj6WicUOr5krYAb2dQE3NVRFwx9v91wJ8ApwBfBX4uImpbXcwlgZqN/9Db\nnLnNK2Np6/kpos5Mr4rPKVkyrSNALeIzrqok8E93Pi3Xvj/01DuOez9Jq4BbgZ9gsGzkHmBbRNyc\n2OcvgF0R8V5JzwdeFREvL5PuSdwmMEdZbQKT9h3dH6+zzVuP2+QrwqQ+ZPxQf11/VQFglM46Ppey\n05ks2nKsynXLsAk4EBGHIuIosBPYOrbPM4FPA0TE36f8v1KuDpqz8R9Y1pw741di4935JrUNtPEH\n1qbqgCrkPc5JJcm8o8GLmtYRYZo8x9bmz3ml3LXzGuBw4vESg8CQtBd4KfBHkl4CnCTp5Ig4UuaN\nszgIzNGknj7jmfh4Zp+3KilPe0OZH2Bbg0wTlD1vaZl+WvfWpsj6vrZ9IGPWOIG9136bG6799rSn\npz15vE7+DcCVkl4BfAb4CvBQsVTm1/o2geV7zir15VlU/WSeH2xWuqaNEyg7dcQkdWQ0bWsYLCPv\ncaYF7vHnNi3TH6kifXW0Q5zw5AOVtAnsvv3sXPte+PSb09oEnge8OSK2DB+/kcEi81ekvYakxwL7\nI2JdmXRP0vogULRhuMgXs+4MadaMetKVfdHjKxoEq854mtRYPq2qpa5G3LxXyy6FlfsMqmoYvub2\nc3Ltu+Xp+9OCwAnALQwahu8BPg9sj4j9iX2eBDwYESHpLcBDEfHmMumepFdBYB79pSfV68+SlmkZ\nRJHxAdPes+jxVW2RgaBI/Xpd3TlnDeZ9UvR7ndy/qiCw6+Czcu3702felPp+wy6if8DDXUR/W9Ll\nwJ6I+IiklwJvZbB0wWeA1w4bkWvRiyDQ1B9U1qCxSVfGaXXAac+ZNogsLS1FVHlOF5X5V9WYXrQk\nlfX5pb1mWsOwzaaq6qAPH8z3eb/4zBtaMVis9UFg+Z6zgOJ1rE0wnmlPCgrTtk3aN/keyW1lG+jq\n6jkyb1V9N7JKflW+fhO/x21RVRD40G3n59r3pRv2tiIIdKZ30LR65Sb+gNL69E/qQTRp26TtSWUz\n4aadw7LqPp6qA4wt3nLHhld1piQwbpYr4C6ZtfdJnobaeZYA6mo0rqp9qI/frbaqqiSw88Bzc+27\n7awvuCQwD1l1p2n68AOdlnlnlYiSme20TLfOBtRJJZ95XA0XbcB35t8/XSsJdOZoimRgXTYaaDZL\nlVGV3WerPs9Vvl7e48xzrpz5N9f4YMuqHI0Tct3aovUlgaSsD7rrGT+kV0+MX0FX2VhZV+ZXd9/8\nPN0JiwYJB4JmqutzmTAvUCt162gSyjSetk0yoxy/nzcwFs0A6wyso3QvIgCYTbOCct3aolMlgaQ+\nDLyZlMFPG6A2nsnOUu8+qatqk3T18++zRQ3sg+6VBDobBPouT++grBJE0fdpcia7qIzCqpM2nibv\n/nVww3ALzVq10OTMoYq0TWpALpqWpp6rqqfCaOpx9kVVAx3L6Noaw70qCUy7ak2rRpllDp42K9Mv\nvysZ5KTPvMuffdNM+hwW+V3rWkmgV0Fgkjq6ktnizZppO7OvX9Y0G6NtdQ0ULKtN3T/zaP2I4TJr\nDJftirjojKLOSd+a+OMbV3SEc1ZVQp8HFy5CmdHjZZ5f1SyiV+zbkmvfy869phUjhrtVrpnBqD57\n1h/8IjPLvmdSyc9tWjVfkUbwvp/XOs16sTVtPq15Wka5blkkbZF0s6RbJV2W8v+1kj4t6TpJeyW9\nqM7j6W0QqKJ3TJ6r0LpNGyGc1Iar+zrlOUeu+7dpVmJVrlsaSauAK4HNwLnAdknjS5X9BvCBiLgA\n2A68o8bD6V8QqLquMRkIFp3JTsu4imRsbckEy1TjNeEz65OqqlwX/Zktx6pctwybgAMRcWi4UMxO\nYOvYPivA44f3n8hgjeHa9KZhuM6MetFfyqS8aw8Uea0mHd/IeHtOkeDexOPpsq6d75KjgdcAhxOP\nlxgEhqTLgd2Sfgl4DPCCMm84TW+CQNe+iHl09ZjLNOhPG0lt1Rufz6rMqPQmfKezrvLv3PMAh/b8\ny7Snp0WQ8d4524E/jYjfHy5M/14GVUe1WHgQkHQn8HUGRaCjEbFJ0snAB4D1wJ3AyyLi6wtLZEsl\nSz+zrgGQbG9owg8Qjs8Yxu9Py+Sd6S/GeHtLU75PRWV1EV3z3NNZ89zTjz3+v++8JW23JWBd4vEZ\nwN1j+7yaQZsBEXGtpO+V9H0RMTXCzKIJbQIrwPMj4vsjYlQseiPwyYjYCHwaeNPCUpdTUzOWMulq\n6jGN85V9PzQlaJQcMbwH2CBpvaQTgW3ArrF9DjGsApJ0DvBv6goA0IwgII5Px1bg6uH9q4EXzzVF\nM/Bgs/kZb9+Z1O/fAaEbmvS7WmFVrluaiFgGLgV2A/uAnRGxX9Llki4a7varwC9I2gu8D7ikzuNZ\neHUQg/qwj0sK4H9HxLuB0yLiPoCIuFfSKQtN4QyaNLFaFUspNuVYYHKGUHRKbJu/ot+nJgUAgOWS\n478i4hpg49i2HYn7+4EfLvUmBTQhCPxgIqPfLekWjm8oyXQw9h27fzKnsFqn1pDEYpqWaRbRxIa4\nadqQRjtenkbiMt25H4z7OcIDMz13kjZNDpfHwoNARNw7/PuApA8z6C51n6TTIuI+SacD92c9/0zV\n1mheSrInRBsCwqSBc9MaXpukDefaBuruzrtap7Kahy8K74j9M79WUtZAsLZa6NFIeoykk4b3Hwtc\nCNzEoKHkFcPdLgH+ZiEJLKlNASDPyOOmBwCzeSg7bUTTLLokcBrw18P2gO8B3hcRuyV9AfgLSa8C\n7gJ+ZpGJ7LpJ1T9pQaGJwaBNwbbPmva9mcVDK92aRXShQSAi7gDOT9n+IDWPkrPpsjL7Jv6QF5XB\n5lmjosj+1nxtWj84j25VbjVQEzPMSapabWwR6j7XaRm65Vd2wsamWA7lurWFg0DN2pahTlroo0uq\nWD1tWsDs8vkrq22/i6Qys4g20aLbBKxBklUVfcjAilTN5O09laXNmV4VuvR96loX0faEK6tNcvRt\nV6dXzjquKo41z8pkXTyneU06P20Mjiso160tXBLoubZmTrMuopM27XTRxt2sfdrSk2oeiiz32TYu\nCZi1lOvvm6Htpc2HVk7IdWsLlwR6qs0/wjwmVf2Mzy9UZOqCvIvUT0pTG6tAiuj6d6tNVT15uCRg\nnZOnjn7SfmkzwmbVX+cdaZ2Vhraadu4mafs5KDmVdOM4CFgnzdLrJ/ncKq9mu35l3Ldp1MsGAUlb\nJN0s6VZJl6X8/22Srpd0naRbJD1Y5/E4CNSsiT+MJqaparP0BJrH6Oi2n/tJPcn60lW2TBCQtAq4\nksHKYecC2yWdndwnIl4/XGTrAuCPgL+q83gcBKxTZllzePz5yZ4+Xci0rFolSwKbgAMRcSgijgI7\nGSyilWU78OcVH8IjOAj0TNuvRCepYtH58e1VnK8uBZJkcEwbB5F1rF06Bw/Fqly3DGuAw4nHS8Nt\nx5G0DngqgyV2a+PeQXPgScOqN2l6i7xdQfM2bqa9ftHPtEuf/6xtJm0fHzBSstE37clZi2htAz4Y\nEbkX2ZqFg4B1yrTMtopG37ZnYovW9gF0WUHgyN7DHNm7NO3pS8C6xOMzgLsz9t0GvKZo+opyEJiT\nJpQG2vzDm2Te57XtmVidshqKkyvttV1WEHjCeet4wnkP5+93vufatN32ABskrQfuYZDRbx/fSdJG\n4IkRkfoiVXKbwBx14QfQdnXOFdR303oKdeW8RSjXLf25sQxcCuwG9gE7I2K/pMslXZTYdRuDRuPa\nuSQwZ4tc/KRLZj2HXbkabZrxUdVduvIfV3bEcERcA2wc27Zj7PHlpd6kAJcEFmDeP4wu/hDLmDTK\nt60zW87bpNHXXR885hHD1ipd/SHWLU8Ds6Xr+rlZXlmV69YWrg5aoL5MKFa1Ks/XpG6L/lxm0/UL\nj6z6/rZyEFiQqupMs17HGVgxdZ6vJvQMs+q0qaonj/aUWTpq1sxhfA6XtP9bffKeX38O3ROR79YW\nDgILNEsDWtEFObp2BdqE4ylyZd+VUbL2sK4tL+kgYK3RhAAA7mYK/Q5qZcYJNJGDwILV2Z+6zz/U\nJkkrDfizaa+udRF1w7C1QlNKAUmzNPi6R1j7ray0J4PPwyWBBvHVYbomZphVpGnUvjNeQvD3oNlc\nHWSVS15ReiWrdsl7fqetM5wWDJpofNGdPupadZCDgNmcjE9JUXdG2tRA0nbuImq1qPoqcNr8OG3R\nhrSW+bymLYBT9rWrCgSupnqYq4OsFZL1zVlLAVo58whQZaoJq5o224MRH6lsEJC0RdLNkm6VdFnG\nPi+TtE/STZLeW9vB4CDQWLOsl5uVKblHSnvNmvmXyaS7PgtoWZHzlkbSKuBKYDNwLrBd0tlj+2wA\nLgP+XUQ8C3hdDYdxjINAQ1W9fm0yELQpGHQ9I8pzfPNcDtMBYLpYUa5bhk3AgYg4FBFHGSwcs3Vs\nn18A/mdEfAMgIv6ltoPBQaB1ymTgk6qImqrpaVx0dU2VrzltZTAbKFkdtAY4nHi8NNyW9Axgo6R/\nlPRZSZtrOIxjHAQaqu4rsbaUDHxFmt+kc9Xn+aaqVrJ3UFp0GN/7e4ANwI8CPwu8W9LjKzuAlDez\nhsqqy59W91/09cdf2xnvbIr0n2/iOR599nnS1ufvSdZV/nf23c53vnzHtKcvAesSj88A7k7Z558j\nYgW4U9ItwFnAF2dK8BQuCdgj9PWHXUYTrpwX8bk14bgXIpR6e/Qzz2T1xS84dsuwB9ggab2kExks\nKL9rbJ8PAz8OIOn7GASA22s6GgeBNlhEH+2m/MCbko68pn1OdXyWszT8Zilyvvt6wVCmOigiloFL\ngd3APmBnROyXdLmki4b7fBz4qqR9wKeAX42II3Udj6uDWmTa6mFVjgXo6w+8rGmDv8pIq26qMgCM\n9LmqJ5eSo4Ej4hpg49i2HWOPfwX4lXLvlI9LAi3mH2pzTOol1LTRtlUNIuurkl1EG8clgZbzQLDm\nyJoYrqv6WmJo05QQebgk0CF9/EE2TXKw1Tw+j7LjRsq8b2+/b2WGDDeQg0BHzPKDbEPpoW0Zzbwz\nx0W1/bRlnEk9lPPWDg4CHVO08a/pP+Impy2vKs/xPAKMG5Cn6FhJwG0CPTPP+eyr0NYFTNqY5ln0\nMhC0KIPPQ9Gm1Q/GSIoX6OJFJ6ORsn6c491I6+zSWLW2Z6xVns9kSa6q1+16hv7J+CBRslVXUqy7\n6opc+9716stKv988uDqow0aZZtOrfPqgDZlr3ukiJj3uhQZWB0l6Tsq2n8rzXAeBjkrOFlp0DqK8\n/7f2mvWzzSpd9uq7kjFtxHG3+XqXpGeNHkjaDvxGnic6CPTErFeiTfpxt+FqOk3dU0bP8hlVtYRp\nWz+TMhT5bnN2MXC1pHMk/QLwGuDCPE90EOi4Pv5ILZ+iwaPMlBWd0sDqoIi4ncFkdB9iEBAujIiv\n53muewfZVE1anzhPb6E+LY5S9vMoejWftv+ivxNz16C2Xkk38ciQsxo4AficJCLi2dNew0HACimz\nklZVkoEgz6ydiwwETcsg085bnkDQhM+9MZrVofKisi/gIGCtNMvsmV0vFdTF523MyqIT8LCIOAQg\n6UxgKSK+K+n5wLOB9+R5jVxtApKO6xibtq1qkrZIulnSrZIuq/v9LL82Zgy+in1kQPT5mFHJ3kHT\n8jVJl0i6X9J1w9urcqTqQ8CypA3AHwNrgffnOZy8DcMvTNn2opzPnYmkVcCVwGbgXGC7pLPrfE+z\nKnUhk21jsK9bmd5BBfK1nRFxwfD2JzmStRIRDwEvAf4oIt4APDnP8UwMApL+y7DhYaOkGxO3O4Ab\n87xBCZuAAxFxKCKOAjuBrTW/Z2f5xzxfTQ4AeUsDZXsPdVa53kF587Wirc9Hh2MDXg58ZLjtUXme\nOK1N4P3Ax4C3Am9MbP9mRDxYMJFFrQEOJx4vMTiBVlBvfpwN0eQAMC4tEOT5vjSpx1jL5M3XXiLp\nR4BbgddHxNKU130l8J+B34qIOyQ9DXhvngRNDALDfqZfB7YDSDoV+F7gJEknRcRded5kRmmR8Lj4\nejD2Hbt/MqewWqfWmCRLcr1yN4x/jkW74TbVg3E/R3ig8tfNqur5zoHb+NfbDk59esq28VfcBbw/\nIo5K+kXgauAnJr1oRHwZ+KXE4zuA356WGMjZO2g4B8XbgKcA9wPrgf0M6rTqsgSsSzw+A7h7fKcz\nVWcSbJo2dR2cR4moDecBqus6O22iwqonuStitU5lNQ9fFN4R+6t54YxG30dvOItHbzjr2OOvXfOJ\ntN2m5mtji8q/C5jaCUfSWQxqbJ7J4EJ99FpPn/bcvA3DbwGeB9waEU9jEJWuzfncWe0BNkhaL+lE\nBqPhdtX8np1Td8bn7pf9lVaCSNOWwJjbSs5buqn5mqTTEw+3Al/Okao/Bf4X8BDw7xl0D81VHZQ3\nCByNiK8CqyStioi/A56b87kziYhl4FJgN7CPQWt5RaHcqtKm9QnqzozaltkVHVk9vprYpONNa2do\n+vcjrzK9g7LyNUmXSxoN/PolSV+SdP1w31fkSNajI+JTDJYHOBQRbwb+Q57jyTtY7GuSTgI+A7xP\n0v3At3I+d2YRcQ2wse736bJ5j5jt6/wyXTvOaSW8tBlq00ZyT5q5NuucNb6tqeSI4bR8LSJ2JO7/\nGvBrBV/2X4fdTw9IuhT4CnBSnifmLQlsBb4D/DJwDXAQyDVXtS1WV66+qtLozGUB8ky7kWeuoNG2\ntP2zTAowjdbACeSA1wGPYdA4/Bzg5xh0F50qVxCIiG9FxHJEPBQRV0fEHw6rh8wmamIQanwmM2d5\nzsekjH/Sc6YFhbTvRxO/M0kNnUo6gD9j0L7wXOAZDBqVp5pYHSTpm6THNAEREY8vlk7ro8YX7612\nnepA0KBZRBPeB7wBuImCsxtNLAlExOMi4vEpt8c5AFgRnfjxd9ioPn8evckmDU5rxQVDM6uDHoiI\nXRFxx7Bh+NBocrlpPItoxy16KmVrj6xAMM9MOWvMQZMCgxo0i2jCDknvBj4FfHe0MSL+atoTvbKY\nmR0z7Uo9S5VrFs+yzzw1tE3glcD5wBYGnXZ+ipxrDbgk0AMuDRzP52Syad03y1ydj3cjXeSo4pk0\na1GZkR+IiJm607skYHPRtB+4A8Bskhn3SNEMPKvHUGs+k2a2CXxW0jNneaKDQE8sejnIpqk6Ta3J\nwGpQxblsU++hhlYHPQ/YK+mW4XT/N0nKNd2/q4PMbO68YH3ltsz6RJcEzKxSo6v5NlzVz6SB1UHJ\nbqFFu4g6CPTIIq60mnp119kMqgGSje5dmjhuRCv5bm3hIGBWs65lgtN0/ngbWBIow0HArEadzxCH\nph1nl85D2YZhSVsk3SzpVkmXTdjvYkkrki6o4zhGHAR6pqnVM9Z+k6p+srqElg0OCwkuJUoCw+me\nrwQ2M1iZcbuks1P2Own4r9S/eJeDgFld2rT0ZhmLvMpfxLktWRLYBBwYNtweBXYymKp/3P9gsKzk\nd1P+VykHgR4qMud7F9WRac0y1UKfLLKBuPL3LdcmsAY4nHi8NNx2jKTzgTMi4qNVJjuLxwn0mKdO\nKK/VI19LmuU4F3Fuqr7gKdnzJ20e6mMhQ5KA3wcumfKcyjgIWG2aWtpw8GuWVs0bBJlX+d+66za+\ndddt0569BKxLPD4DuDvx+HEM2gr+fhgQTgf+RtJPR8R1syZ5EgcBsxllZVwOMh2XEQQeu3YDj127\n4djjf/mn3Wm77QE2SFoP3ANsA7Yfe+mIbwCnjh5L+jvg9RFxfQUpT+U2gZ5r1RVYg/T5vNVRv9+m\noFmmYTgiloFLgd3APmBnROyXdLmktKmfA1cHWRs1PZMse7U+qQqjTRlaEV09rsJKDgSLiGuAjWPb\ndmTs++Pl3m06lwTMatT0YJjXPAJAW4JMQ2cRnZmDgJm1VtsGizWRq4N6ruofUVeufO1h88xoZ1mg\nZt7adJWfh4OAVcKZv1Wl8V1GOxYEXB1kpTX6BztBW9PdB01uH+ham4BLAmYFTQseDi7lNfoctiiD\nz8NBwEpp9I/VrA4dCwKuDrJSmlxsN6tD16qDHASs11ySmW4R56jRFxfuImpdUcUPrW+ZaN+O147X\npvWD83AQsJk5Q7Q+alNVTx6uDrKZdCkAdOlYuqLRn0nHqoMcBMxyanTGVLM+H/txOhYEXB1kloMz\nwcn6NKNq2eogSVuAtzO4CL8qIq4Y+/8vAq8FloFvAv8pIm4u967ZXBIwMyuiRElA0irgSmAzgxXE\ntks6e2y390XEsyPi+4HfZbDcZG0cBKywPl0Vb37Keb063iyTrugnnZ9Zzl3Tz7cict0ybAIORMSh\niDgK7AS2JneIiP+XeHgSUGt/JFcH2VSjH2XjJ/YqYbTITFePr6zkIjxlMvYuVA+V7CK6BjiceLzE\nIDA88j2k1wCvBx4F1LqwjEsCNlWZH3+bdP34ypq1VDTv81p7oCnXMJy2VORxe0fEOyJiA3AZ8N9L\np3kClwR6qgtXZNY9VQSMuoNOVsPwN+69jW/ee3Da05eAdYnHZwB3T9j/A8A7CySvMJcELBcHDavC\nvEoFH7/7hvq+sxlX/o8/bQNrztt87JZhD7BB0npJJwLbgF3JHSRtSDy8CLi16kNIchCwXFxVYm1S\n5/e1zARyEbEMXArsBvYBOyNiv6TLJV003O1SSV+SdB3wOuCS2g4GVwfZFM78rWrz+E6NSgG1NPaX\nHCcQEdcAG8e27Ujcf125dyjGQcBSOfO3Lqjje+y5g8zMGq7W6qCVyHVrCwcBM1uYVnY48NxB1gXJ\nwT/j283mpewgtHHJ16rru+z1BKwTWnkFZtYELbrKz8NBoIccAKxpqi4F1MkNw9Z6rvKxJmlTAAAg\nIt+tJVwS6KGskoCDg7XRvEu2bhOwTnIAsLab13fY1UEVkbRD0pKk64a3LYn/vUnSAUn7JV24qDT2\nidsJrM3mehHTseqgRbcJvC0iLhjergGQdA7wMuAc4EXAOySlTb9qFXMgsDaadym2zNxBTbToIJCW\nuW9lMKnSQxFxJ3CAlEUXrHquEjLLoWODxRYdBF4raa+kd0t6wnDb+Mo7Xxlusxo5AJjl07WSQK0N\nw5I+AZyW3MQgRv468A7gNyMiJL0F+D3g58m58s7Iwdh37P7JnMJqnVpByrtrNKuiq36s6x6M+znC\nA9W/cIvmBcqj1iAQES/Mueu7gL8d3l8C1ib+N3HlnTN17myJ66m0AOBSgHXRap3Kah6+KLwj9lfy\numW7iA47wbydQU3MVRFxxdj/f5nBBfFR4AHgVRFx+LgXqsgiewednnj4EuBLw/u7gG2STpT0NGAD\n8Pl5p8/MLFWJ3kGSVgFXApuBc4Htks4e2+064DkRcT7wIeB3azyahbYJ/I6kGyXtBX4M+GWAiPgy\n8BfAl4GPAq+JaFF/qxYYv/J31ZBZfiXbBDYBByLiUEQcBXYy6AxzTET8Q0T86/DhtdTcJrqwwWIR\n8fIJ/3sr8NY5Jqd33C5gNqNyl6TjHV+WmNz78dXAx0q94xQeMdxTyblWHAzM8lNGxcSRB2/na0du\nn/r0lG2L2TtDAAAIeElEQVSpLyjp54DnMKgpqY2DQE+NqoSSa7Emt1eplnVezRYlo2H45Cc+nZOf\n+PRjj++841Npuy0B6xKPUzu+SHoB8CbgR4fVRrVZ9DgBW7Bk5lxXRu0AYF2iiFy3DHuADZLWSzoR\n2MagM8zDry99P/BO4Kcj4qu1HgwuCRjOpM0KKTFOICKWJV0K7ObhLqL7JV0O7ImIjwC/AzwW+Mvh\nlDmHIuLFFaQ8lYOAmVkBZUcDD+dJ2zi2bUfift7xVZVwEDAzK6JjPdYdBMzMCvCiMmZmfeaSgJlZ\nj3UrBjgImJkVMaH7Zys5CJiZFbHsIGBm1lsuCZiZ9ZmDgJlZjzkImJn1mMcJmJn1l9sEzMz6zEHA\nzKzHVrpVH+T1BMzMiljJecsgaYukmyXdKumylP//iKQvSjoq6SU1HMEjOAiYmRVQZlEZSauAK4HN\nwLnAdklnj+12CLgEeF+dxzHi6iAzsyLKtQlsAg5ExCEASTuBrcDND7983DX831waH1wSMDMrYiXy\n3dKtAQ4nHi8Nty2MSwJmZkVklAS++u27ePA7d017ttJesWySynAQMDMrIiMIPOnRa3nSo9cee3zw\nwc+m7bYErEs8PgO4u8LUFeYgYGZWxHKpLqJ7gA2S1gP3ANuA7RP2Tys5VMptAmZmRcRKvlvaUyOW\ngUuB3cA+YGdE7Jd0uaSLACQ9V9Jh4GLgnZJuqvNwXBIwMyui5IjhiLgG2Di2bUfi/heAtePPq4uD\ngJlZEdk9f1rJQcDMrAjPHWRm1mMOAmZmPba8vOgUVMpBwMysCJcEzMx6zEHAzKzH3DvIzKy/ImMg\nWFs5CJiZFeGSgJlZj7lNwMysx9xF1Mysv6JjC807CJiZFdGx6iBPJW1mVkS55SWRtEXSzZJulXRZ\nyv9PlLRT0gFJ/yxpXdrrVMVBIMODcf+ik3Acpykfpym/JqariWl6hBLrCUhaBVwJbAbOBbZLOnts\nt1cDD0bEWcDbgd+p8WgcBLIc4YFFJ+E4TlM+TlN+TUxXE9OUFCuR65ZhE3AgIg5FxFFgJ7B1bJ+t\nwNXD+x8EfqKWAxlyEDAzK6JESQBYAxxOPF4abkvdZ7gS2dckra76MEbcMGxmVkCU6yKatmbweLFh\nfB+l7FMZRYtbuiW1N/FmNncRUWrhdkl3Autz7n5fRJw+9vznAW+OiC3Dx28cJCuuSOzzseE+n5N0\nAnBPRJxaJt2TtLokUPYDNTMrIiKeWvIl9gAbJK0H7gG2AdvH9vlb4BLgc8DPAJ8u+Z4TtToImJm1\nSUQsS7oU2M2gTfaqiNgv6XJgT0R8BLgK+DNJB4CvMggUtWl1dZCZmZXT+95Bki6W9CVJy5IuSGxf\nL+nbkq4b3t6R+N8Fkm4cDvZ4+7zSNPzfm4aDSPZLujCxfeIAlBrSuEPSUuL8bJmWxnmY93mYkI47\nJd0g6XpJnx9uO1nSbkm3SPq4pCfUnIarJN0n6cbEtsw0SPrD4ee2V9L5c0xTI79LvRERvb4BG4Gz\nGNS7XZDYvh64MeM5nwM2De9/FNg8pzSdA1zPoBrvqcBtDHoOrBreXw88CtgLnF3zedsBvD5le2oa\n5/RZzv08TEjL7cDJY9uuAP7b8P5lwG/XnIYfBs5Pfo+z0gC8CPg/w/v/Frh2jmlq3HepT7felwQi\n4paIOEB6163jtkk6HXhcRHx+uOk9wIvnlKatwM6IeCgi7gQOMBh8kmcASh3SzllWGudhUechzSg4\nJyUHAV1Nxd+bcRHxj8CRKWnYmtj+nuHzPgc8QdJpc0oTNO+71Bu9DwJTPFXSFyX9naQfHm5bw2CA\nx0jaYI+6jA80+cpwW54BKHV47bDq4N2JaoWsNM7Dos5DmgA+LmmPpJ8fbjstIu4DiIh7gVMWkK5T\nx9Iw6nq4yM8Nmvdd6o1e9A6S9AkgeVUzGnzx6xHxtxlPuxtYFxFHhvXyH5b0TPIN9qgrTVnvnRbM\nS7f4T0oj8A7gNyMiJL0F+D3g5yekcR4W+d7jfjAi7pV0CrBb0i0LTEseizx3Tfwu9UYvgkBEvHCG\n5xxlWGyNiOskHQSeweDqcm1i1zMYBIza0zThvQWsS9leSoE0votB3+ZJaZyHJWo4D7MYXmUTEQ9I\n+jCDaoz7JJ0WEfcNqxUXMVNaVhoW9rlFRHKyoKZ8l3rD1UGPdOzKQ9L3DWf8Q9LTgQ3A7cMf9zck\nbZIk4OXA38wjTcAuYNtwqtmnDdP0eRIDUCSdyKBf8a4a0zRqGxl5CfClKWmch7mfhzSSHiPppOH9\nxwIXAjcN0/KK4W6XUO/35lhyOP47NErDKxJp2MXguzwa1fq1UbVR3Wlq6HepPxbdMr3oG4PGucPA\ndxiM4PvYcPvoy3g98AXgJxPPeQ6DH/UB4A/mlabh/97EoJfEfuDCxPYtwC3DNL1xDuftPcCNDHrg\nfJhBfffENM7p85zrechIw9OG5+X64ffkjcPtq4FPDtP3CeCJNafj/QyunL8L3AW8Ejg5Kw0Mpji+\nDbiBRK+0OaSpkd+lvtw8WMzMrMdcHWRm1mMOAmZmPeYgYGbWYw4CZmY95iBgZtZjDgJmZj3mIGCt\nJOmbi06DWRc4CFhbeYCLWQUcBKz1JP2upJuGi7i8bLjtx4azv/7lcEGSP1t0Os2aqBcTyFl3SXop\n8OyIeJakU4E9kv5h+O/zgWcC9wL/JOkHI+Kzi0qrWRO5JGBt90PAnwNExP3A3wM/MPzf5yPinhjM\njbKXwepUZpbgIGBtNz7nfPLxdxP3l3HJ1+w4DgLWVqPM/jPAf5S0ariAy4/g6YbNcvOVkbVVAETE\nXw/nv78BWAHeEBH3SzonbX8zeyRPJW1m1mOuDjIz6zEHATOzHnMQMDPrMQcBM7MecxAwM+sxBwEz\nsx5zEDAz6zEHATOzHvv/WRYGXX8hxJ4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118006240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Extract the area variable\n",
    "area = xr.open_dataset('area.nc').load()['cell_area']\n",
    "dom_ds['area'] = area\n",
    "\n",
    "# write the domain file\n",
    "dom_ds.to_netcdf(dom_file)\n",
    "dom_ds.close()\n",
    "\n",
    "# document and plot the domain\n",
    "print(dom_ds)\n",
    "dom_ds.mask.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!rm area.nc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maurer 1/8 deg. Parameters\n",
    "\n",
    "Finally, we'll repeat the same steps for the Maurer 1/8 deg. parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading ./Maurer_0.125_NLDAS/soil/us_all.soil.wsne\n",
      "reading ./Maurer_0.125_NLDAS/snow/us_all.snowbands.wsne\n",
      "reading ./Maurer_0.125_NLDAS/veg/us_all.veg.wsne\n",
      "reading ./Livneh_0.0625_NLDAS/LDAS_veg_lib\n",
      "Calculating grid size now...\n",
      "found 462 unique lons\n",
      "found 222 unique lats\n",
      "creating 2d coordinate arrays\n",
      "Created a target grid based on the lats and lon in the soil parameter file\n",
      "Grid Size: (222, 462)\n",
      "gridding params now...\n",
      "creating 2d coordinate arrays\n",
      "48547 masked values\n",
      "writing var: run_cell\n",
      "writing var: gridcell\n",
      "writing var: lats\n",
      "writing var: lons\n",
      "writing var: infilt\n",
      "writing var: Ds\n",
      "writing var: Dsmax\n",
      "writing var: Ws\n",
      "writing var: c\n",
      "writing var: expt\n",
      "writing var: Ksat\n",
      "writing var: phi_s\n",
      "writing var: init_moist\n",
      "writing var: elev\n",
      "writing var: depth\n",
      "writing var: avg_T\n",
      "writing var: dp\n",
      "writing var: bubble\n",
      "writing var: quartz\n",
      "writing var: bulk_density\n",
      "writing var: soil_density\n",
      "writing var: off_gmt\n",
      "writing var: Wcr_FRACT\n",
      "writing var: Wpwp_FRACT\n",
      "writing var: rough\n",
      "writing var: snow_rough\n",
      "writing var: annual_prec\n",
      "writing var: resid_moist\n",
      "writing var: fs_active\n",
      "writing var: cellnum\n",
      "writing var: AreaFract\n",
      "writing var: elevation\n",
      "writing var: Pfactor\n",
      "writing var: Nveg (222, 462)\n",
      "writing var: Cv (12, 222, 462)\n",
      "writing var: root_depth (12, 2, 222, 462)\n",
      "writing var: root_fract (12, 2, 222, 462)\n",
      "writing var: LAI (12, 12, 222, 462)\n",
      "writing var: overstory (12, 222, 462)\n",
      "writing var: rarc (12, 222, 462)\n",
      "writing var: rmin (12, 222, 462)\n",
      "writing var: wind_h (12, 222, 462)\n",
      "writing var: RGL (12, 222, 462)\n",
      "writing var: rad_atten (12, 222, 462)\n",
      "writing var: wind_atten (12, 222, 462)\n",
      "writing var: trunk_ratio (12, 222, 462)\n",
      "writing var: albedo (12, 12, 222, 462)\n",
      "writing var: veg_rough (12, 12, 222, 462)\n",
      "writing var: displacement (12, 12, 222, 462)\n"
     ]
    }
   ],
   "source": [
    "soil_file = os.path.join(mdpath, 'soil', 'us_all.soil.wsne')\n",
    "snow_file = os.path.join(mdpath, 'snow', 'us_all.snowbands.wsne')\n",
    "veg_file = os.path.join(mdpath, 'veg', 'us_all.veg.wsne')\n",
    "vegl_file = os.path.join(ldpath, 'LDAS_veg_lib')  # from livneh\n",
    "out_file = os.path.join(mdpath, 'nldas_0.125deg.vic_5.0.0_parameters.nc')\n",
    "\n",
    "cols = gp.Cols(nlayers=3,\n",
    "               snow_bands=5,\n",
    "               organic_fract=False,\n",
    "               spatial_frost=False,\n",
    "               spatial_snow=False,\n",
    "               july_tavg_supplied=False,\n",
    "               veglib_fcan=False,\n",
    "               veglib_photo=False)\n",
    "n_veg_classes = 11\n",
    "root_zones = 2\n",
    "vegparam_lai = True\n",
    "lai_src = 'FROM_VEGPARAM'\n",
    "\n",
    "\n",
    "# ----------------------------------------------------------------- #\n",
    "\n",
    "# Read the soil parameters\n",
    "soil_dict = gp.soil(soil_file, c=cols)\n",
    "\n",
    "# Read the snow parameters\n",
    "snow_dict = gp.snow(snow_file, soil_dict, c=cols)\n",
    "\n",
    "# Read the veg parameter file\n",
    "veg_dict = gp.veg(veg_file, soil_dict,\n",
    "                  vegparam_lai=vegparam_lai, lai_src=lai_src, \n",
    "                  veg_classes=n_veg_classes, max_roots=root_zones)\n",
    "\n",
    "# Read the veg library file\n",
    "veg_lib, lib_bare_idx = gp.veg_class(vegl_file, c=cols)\n",
    "\n",
    "# Determine the grid shape\n",
    "target_grid, target_attrs = gp.calc_grid(soil_dict['lats'], soil_dict['lons'])\n",
    "\n",
    "# Grid all the parameters\n",
    "grid_dict = gp.grid_params(soil_dict, target_grid, version_in='4', \n",
    "                           vegparam_lai=vegparam_lai, lai_src=lai_src,\n",
    "                           lib_bare_idx=lib_bare_idx, \n",
    "                           veg_dict=veg_dict, veglib_dict=veg_lib, \n",
    "                           snow_dict=snow_dict, lake_dict=None)\n",
    "\n",
    "# Write a netCDF file with all the parameters\n",
    "gp.write_netcdf(out_file,\n",
    "                target_attrs,\n",
    "                target_grid=target_grid,\n",
    "                vegparam_lai=vegparam_lai,\n",
    "                lai_src=lai_src,\n",
    "                soil_grid=grid_dict['soil_dict'],\n",
    "                snow_grid=grid_dict['snow_dict'],\n",
    "                veg_grid=grid_dict['veg_dict'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1/8 deg CONUS domain file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dom_ds = xr.Dataset()\n",
    "\n",
    "# Set global attributes\n",
    "dom_ds.attrs['title'] = 'VIC domain data'\n",
    "dom_ds.attrs['Conventions'] = 'CF-1.6'\n",
    "dom_ds.attrs['history'] = 'created by %s, %s' % (user, now)\n",
    "dom_ds.attrs['user_comment'] = 'VIC domain data'\n",
    "dom_ds.attrs['source'] = 'generated from VIC CONUS 1.8 deg model parameters, see Maurer et al. (2002) for more information'\n",
    "\n",
    "dom_file = os.path.join(opath, 'domain.vic.conus0.0125deg.%s.nc' % now.strftime('%Y%m%d'))\n",
    "\n",
    "# Get the mask variable\n",
    "dom_ds['mask'] = xr.DataArray(target_grid['mask'], coords={'lat': target_grid['yc'],\n",
    "                                                                 'lon': target_grid['xc']},\n",
    "                                    dims=('lat', 'lon', ))\n",
    "\n",
    "# For now, the frac variable is going to be just like the mask\n",
    "dom_ds['frac'] = dom_ds['mask'].astype(np.float)\n",
    "dom_ds['frac'].attrs['long_name'] = 'fraction of grid cell that is active'\n",
    "dom_ds['frac'].attrs['units'] = '1'\n",
    "\n",
    "# Set variable attributes\n",
    "for k, v in target_attrs.items():\n",
    "    if k == 'xc':\n",
    "        k = 'lon'\n",
    "    elif k == 'yc':\n",
    "        k = 'lat'\n",
    "    dom_ds[k].attrs = v\n",
    "    \n",
    "# Write temporary file for gridarea calculation\n",
    "dom_ds.to_netcdf('temp.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cdo gridarea:   0%\b\b\b\b  0%\b\b\b\b  1%\b\b\b\b  2%\b\b\b\b  3%\b\b\b\b  4%\b\b\b\b  5%\b\b\b\b  6%\b\b\b\b  7%\b\b\b\b  8%\b\b\b\b  9%\b\b\b\b 10%\b\b\b\b 11%\b\b\b\b 12%\b\b\b\b 13%\b\b\b\b 14%\b\b\b\b 15%\b\b\b\b 16%\b\b\b\b 17%\b\b\b\b 18%\b\b\b\b 19%\b\b\b\b 20%\b\b\b\b 21%\b\b\b\b 22%\b\b\b\b 23%\b\b\b\b 24%\b\b\b\b 25%\b\b\b\b 26%\b\b\b\b 27%\b\b\b\b 28%\b\b\b\b 29%\b\b\b\b 30%\b\b\b\b 31%\b\b\b\b 32%\b\b\b\b 33%\b\b\b\b 34%\b\b\b\b 35%\b\b\b\b 36%\b\b\b\b 37%\b\b\b\b 38%\b\b\b\b 39%\b\b\b\b 40%\b\b\b\b 41%\b\b\b\b 42%\b\b\b\b 43%\b\b\b\b 44%\b\b\b\b 45%\b\b\b\b 46%\b\b\b\b 47%\b\b\b\b 48%\b\b\b\b 49%\b\b\b\b 50%\b\b\b\b 51%\b\b\b\b 52%\b\b\b\b 53%\b\b\b\b 54%\b\b\b\b 55%\b\b\b\b 56%\b\b\b\b 57%\b\b\b\b 58%\b\b\b\b 59%\b\b\b\b 60%\b\b\b\b 61%\b\b\b\b 62%\b\b\b\b 63%\b\b\b\b 64%\b\b\b\b 65%\b\b\b\b 66%\b\b\b\b 67%\b\b\b\b 68%\b\b\b\b 69%\b\b\b\b 70%\b\b\b\b 71%\b\b\b\b 72%\b\b\b\b 73%\b\b\b\b 74%\b\b\b\b 75%\b\b\b\b 76%\b\b\b\b 77%\b\b\b\b 78%\b\b\b\b 79%\b\b\b\b 80%\b\b\b\b 81%\b\b\b\b 82%\b\b\b\b 83%\b\b\b\b 84%\b\b\b\b 85%\b\b\b\b 86%\b\b\b\b 87%\b\b\b\b 88%\b\b\b\b 89%\b\b\b\b 90%\b\b\b\b 91%\b\b\b\b 92%\b\b\b\b 93%\b\b\b\b 94%\b\b\b\b 95%\b\b\b\b 96%\b\b\b\b 97%\b\b\b\b 98%\b\b\b\b 99%\b\b\b\b100%\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \b\b \bcdo gridarea: Processed 2 variables ( 0.12s )\n"
     ]
    }
   ],
   "source": [
    "# This step calculates the grid cell area\n",
    "!cdo -O gridarea temp.nc area.nc\n",
    "!rm temp.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:  (lat: 222, lon: 462)\n",
      "Coordinates:\n",
      "  * lat      (lat) float64 25.19 25.31 25.44 25.56 25.69 25.81 25.94 26.06 ...\n",
      "  * lon      (lon) float64 -124.7 -124.6 -124.4 -124.3 -124.2 -124.1 -123.9 ...\n",
      "Data variables:\n",
      "    mask     (lat, lon) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...\n",
      "    frac     (lat, lon) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    area     (lat, lon) float64 1.748e+08 1.748e+08 1.748e+08 1.748e+08 ...\n",
      "Attributes:\n",
      "    title: VIC domain data\n",
      "    Conventions: CF-1.6\n",
      "    history: created by jhamman, 2016-06-08 15:01:23.814110\n",
      "    user_comment: VIC domain data\n",
      "    source: generated from VIC CONUS 1.8 deg model parameters, see Maurer et al. (2002) for more information\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x11d8317f0>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEPCAYAAABIut/fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+8HHV97/HXOwg2iiJR+ZWQgASD5qEgWKRUL1CQBFTw\nItrkXgV/9NbeNvUX+gCsbUxrq7T+oErvxSpSUGxoRSXlYogKYhXBo5AIMYHwK+TITw0UK2jDOZ/7\nx845mWx298ye2dmdmX0/H495ZGd2Zva7J7vv+e53vvMdRQRmZlY/MwZdADMzK4YD3sysphzwZmY1\n5YA3M6spB7yZWU054M3MasoBb2bWI5IukvSQpJ90WOfTkjZJWivpsCLL44A3M+udi4FF7Z6UdBJw\nUEQcDLwTuLDIwjjgzcx6JCK+BzzaYZVTgUuTdW8C9pC0d1HlccCbmfXPbGBLav5nybJCOODNzPpH\nLZYVNl7M04racS9I8kA5ZpZZRLQK0EwO2H/X2Dz6VDebPBQR+3T5MqPA/qn5OcD9Xe4js1IHPMAJ\nOr3n+7wr1nPnAzv/Ry7a79Cev1Yv3RXrOUgLB12Maalq2V3u/ptu2b8VX8n1uptHn2LbAwdlXn/X\nfe9q13YuWtfUAVYBfwJcLuko4LGIeKirgnah9AFflEX7Hco196/bYdnEfNmD3syKMRbjubaX9GXg\nWOC5ku4DlgO7ARER/xgRV0s6WdKdwK+At+UsckdDG/CwPchbBb1D3mz4jOdsDo+I/5FhnWW5XqQL\nQxnwe/L8QRdhWqpabqhu2V3u/htk2bfF2MBeuwhDGfCztNeU66Rr9WWpzWcpd1lVtewud/8Nsux5\na/Bl426SGTQ34ZhZPY0RmacqGMoavJlZK3WrwTvgM3IPG7P6G6vZPardRNMlN9eY1dd4F1MVOOBx\nrdzMGurWBu+AnwbX4s3qaVtkn6rAAT9NDnmz+hlDmacqcMAn3ExjZuORfaoC96LJwT1rzOqlKjXz\nrFyDT3FQmw03N9HYTtweb1YP46HMUxUMdcA7mM0szTX4Lkm6V9I6SbdI+mGybE9JayTdLukaSXsU\nXY52ehXy19y/zgcMs4rbFrtknlqRtFjSRkl3SDq7xfNzJX0rycRrJe1X5PvpRw1+HDg2Il4WEUcm\ny84BvhURC4BrgXP7UI4dFBXGDnmz6spTg5c0A7gAWAQsBJZKOqRptY8D/xQRhwJ/CXysyPfTj4BX\ni9c5FbgkeXwJ8Po+lMPMrKOxmJF5auFIYFNEbI6IbcBKGlmX9mIalVoi4jstnu+pfgR8ANdIGpH0\nB8myvSfuQxgRD0J/R/hvV8t27dtsuI0zI/PUwmxgS2p+NFmWthZ4A4Ck04DdJe1ZxHuB/vSDPzoi\nHpT0fGCNpNuhXAM5ONjNDDr3g1974xOsu/GJTpu32rg56z4AXCDprcB3gZ8BT3VXyuwKD/ikhk5E\nPCLp6zR+xjwkae+IeEjSPsDD7ba/K9ZPPt6T50/rbi/pC5Ic5mb1sDUe5lEe6ek+2zS9APCSV+zO\nS16x++T8pX+/tXmVUWBuan4OcH96hYh4gO01+GcCb4iIX+YqdAeFBrykZwAzIuI/kzdzIrACWAW8\nFTgPOBO4st0+DtLCXGVIB7rD3aw+ZmkvZrG9wndPbMi9z/F83R9HgPmS5gEPAEuApekVJD0X2BoR\nQaNzyRfyvOBUiq7B7w18TVIkr3VZRKyR9CPgXyS9HbgPeGPB5eira+5f1/OrYjsdnHwFrllv/FdM\nPxIjYkzSMmANjfObF0XEBkkrgJGIuAo4FviopHEaTTR/kr/U7RUa8BFxD3BYi+VbgROKfG2oR409\ny3totY5D36x7bU6eZhYRq4EFTcuWpx5fAVyR60W6MNRXstZZHQ5uZv02Fso8VYFHkxyQfoxEmW4q\nahf4rU48V6X2X0RTmA23sZrVeSsb8OlQavUlH3QNtlOA9/PE71T7b/W8h0G2YTXeoRdNFVXy3Qw6\nvPOoUtnLXlYfgKzXxpiReaqCStbgm5sdyhxEZS5bFmUu/0TzkoPeeqXdIGJVVY3DUBv+Yg+3Khzg\nrVpyjkVTOpWswZs1q+qJYiuXnBc6lU7lA97DD1grDnybjqrUzLOqfMA73C0L9wyyLKpy8jSr0ge8\nT6JZL03VvdaGW1XutZpV6QMe2l+K79q7wc69qrJqdSGYQ3+4uQZfEh4C2CD//3/ztp2u+E2v4wNB\n//TzV1fduklWNuAnONyHV79C1uE+POp2JWvlA96GV9EH91ZB7nDvr3a/sIr6f+h0R6csJC0Gzmf7\ncMHnNT2/P437UD8nWefciPhGrhftwAFv1sQhXg6D+HWepwYvaQZwAXA8jTs5jUi6MiI2plb7EHB5\nRHxW0ouAq4EDcxS5Iwe8WcLBXh2tTpDvsm/+/ebsB38ksCkiNgNIWgmcCqQDfhx4dvL4OTTuyVoY\nB7xZDu59U4ystfde1/JzXsk6G9iSmh+lEfppK4A1kt4FPIOCb3zkgDej+4B2f/rBK6IJp1MN/t6R\nR9g88vNOm7c6OkTT/FLg4oj4lKSjgC8B+W483UHpA97dIK1oDvdyGeT3vVM3ydkv34fZL99ncv7f\nL7y9eZVRYG5qfg6Ntvi0dwCLACLiRkm/Jel5EdHxyDFd9eoTZNalPAHtcO+9QVfmxkOZpxZGgPmS\n5knaDVgCrGpaZzNJs0xykvXpRYU7VKAGb1aU5oCeqo+7a+7FGnS4Q76bbkfEmKRlwBq2d5PcIGkF\nMBIRVwHvBz4n6b00Trie2YNit+WAt6HUfPFSJ51ua9hun1ZNeW+mHRGrgQVNy5anHm8AXpnrRbrg\ngLehNdWQwkXVKN3zpvcaf8tNuffjwcbMaiBrrbzbffbrIFE3ef5OvTxQeqgCM2urm6DyTUl6o6cX\nOtXsjk71OlyZWSWV5eD21PgumacqqETAl+U/38yKUZamrHGUeaqCSgS8WV0s2u/QzF0xh0lZKnFj\nocxTFbgN3qxg7YYdbtfv3r1sBscnWc0ss04h7fHmt2u+iGxQv2TcTdLMOsoa0q6pN2S9bWI/VKVt\nPSsHvFkPdRvu6cfDHvRlULcafF8anCTNkHSLpFXJ/D9JujtZdrOkl/ajHGZl0antfZiU7T3XrZtk\nv2rw7wbWs/1OJgGcFRFf69Prm5VKu2Abltp82YJ9Qt2aaAqvwUuaA5wMfL7fr23WT9NpnhlGZX7/\nOYcLLp1+1OA/BXwA2KNp+Uck/TnwbeCciNjWbgdl/kCYpXUacniq4YjrYFDf1bIMNiZpMXA+24cL\nPq/p+U8Cx9FoxXgm8PyImJXrRTsoNOAlvQZ4KCLWSjo29dQ5EfGQpF2BzwFnAx9ptY+7Yj0rPv4L\nAI45eibHHv2MIotsNi1TBXddKilleh/fueEJrr/hSaCRE72QJ+AlzQAuAI6ncSenEUlXRsTkTbcj\n4n2p9ZcBh02/tBnKFNF8y8Ae7lz6G+DNwFPATOBZwFcj4ozUOsfQaI8/pcX2cYJOL9WHyqzZRLhP\ndeHSdPdbFlnfx6D6se+y7yYipp/QkuLk6/808/pXH/OZHV4vucfq8og4KZk/B4jmWnxq/e8DfxER\n355umadSaDt4RHwwIuZGxAto3L7q2og4Q9I+AJIEvB64rchymBWt7pWQPKNk9kOvDoZPxYzMUwuz\ngS2p+dFk2U4kzQUOAK7tScHbGFQ/+MskPY/GXcjXAn80oHKY5Vb3ni9lPXgV8ffO2QbfauN2TSRL\ngK9EkU0o9DHgI+J64Prk8fH9el2zKqrrwWI6msftKfJkdaeAf3TtFh5dO9pp81Fgbmp+Do22+FaW\nAH/cbfm65StZzaytMtbeizz4dQr4PQ6dyx6Hbs/vey+9sXmVEWC+pHnAAzRCfGnzSpIWAM+JiJ12\n0GsOeLMClTEgsypD2SfCvF+/aHKcoyUixpKeMWvY3k1yg6QVwEhEXJWsugRYmbuwGTjgzaahyJ4i\nbp4ZnLxXskbEamBB07LlTfMrcr1IF3w1qVkP9CqU0/u55v51A61FD/pAM4jX95WsZlZ48DaPNjmI\nsKtyd8fpGhuvV53XAW9mpTDocId8bfBl5IA3K5EynNgcZlVpesnKAW9mlij2sqP+K33Au0ZjZVeG\npoVeG9bvXd3Ggy99wJtZb7QK7bLcWaosB0m3wZtZZoMaWbHZVHeQaqcf5S9LuIPb4M1siJTh4NRP\n4+P1Cvh6dfo0q6m8QVumWnKZRSjzVAUOeLOchq2W20tlO/DU7UpWB7yZDUzZDo4R2acqcBu8Wcn1\nopZbtiAtq6o0vWTlGrxZzZU53MvWRJO3DV7SYkkbJd0h6ew267xJ0npJt0r6UpHvxzV4s5zKFlJp\nZQ73MsrT8iJpBnABcDyNOzmNSLoyIjam1pkPnA38TkQ8nty6tDAOeLOcOo32mCdg8xw4qhLsgxop\ns53I103ySGBTRGwGkLQSOBXYmFrnfwH/EBGPA0TEz/O84FTcRGPWA2UK1DKVpWpyNtHMBrak5keT\nZWkvBBZI+p6kGyQtKuitABWowZflSkCzbvX7c+vvSX45e8e0Sv3mPT4NmA/8Nxo36P53SQsnavS9\nVvqANxtG3TRbVD3Yy9RM06kXzZPr7+bJn97TafNRGqE9YQ6NtvjmdX4QEePAvZJuBw4GfjytAk/B\nAW/WY/0K3KoHe9rEexl40HcI+JkvPoiZLz5ocv6xK65rXmUEmC9pHvAAjZtrL21a5+vJskuTE6wH\nA3fnLncbDnizHull4E4VeHUJ94km2IEHeyJPE01EjElaBqyhcX7zoojYIGkFMBIRV0XENZJOlLQe\neAp4f0Q82ouyt+KAHwCfV7Cs0kFfl89Mc5iXJdyBfP0kgYhYDSxoWra8af4s4Kx8r5SNA35A6vSF\nteJV/bNSqhDvIGc3ydJxwA9AmX6SmhWhqp/vug1V4IAfINfirW6qGuyTKjKIWFYO+AFxsFvdpMO9\n+fNdneB3Dd7MbCdZ7/laaq7B95drumbVUrlQT3PAdy8ZZe1HwGhEnCLpAGAlsCdwM/CWiHiqH2Ux\ns96qdKA3qVsvmn4NNvZu4Kep+fOAT0TEAuAx4B19KoeZWXvRxdQnko5osex1WbYtPOAlzQFOBj6f\nWvx7wBXJ40uA/150OcysGLVqRg1ln/rnc5JeMjEjaSnwoSwb9qOJ5lPAB4A9ACQ9F3g0GWwHGoPv\n7NeHcpiZdaRytsGfDnxF0v8EXgmcAZyYZcNCA17Sa4CHImKtpGMnFrNzX6S2f9b5+z6Nt5z1EADH\nHD2TY49+RhFFNbNpGlQb/NZ4mEd5pLc7LWHAR8TdkpbQGKhsC3BiRDyZZduia/C/C5wi6WRgJvAs\n4HxgD0kzklp8qyE1Jx2khSx/v8+/mtmOZmkvZrHX5Pw9sSH/Tkt0JaukW9nxkDML2AW4SRIR8dKp\n9lFowEfEB4EPAkg6BjgrIt4s6XLgjcDlwJnAlZ324ys+zawvylWDf23eHQzqln3nAO+TdAeNo9JF\nAyqHmdl2411MBYuIzcn9XZ8GPJg8PpDGfV7/I8s+MgW8pPOyLOskIq6PiFOSx/dExCsi4oUR8fsR\nsa2bfZlZf7VrZ69TH3ggdy8aSYslbZR0h6SzWzx/pqSHJd2cTG/PUKorgDFJ84F/BPYHvpzl7WSt\nwb+6xbKTMm5rZhU3LE2kiuzTTts2Lui8AFgELASWSjqkxcusjIjDk+kLGYo1nlwIehrwmYj4ALBv\nlvfTMeAl/e+koX+BpJ+kpnuAn2R5ATOrpqluzFG72jvkvdDpSGBT0rSyjcbV+qe2WK/bM7nbkr7v\nZwBXJct2zbLhVDX4LwOvA1Yl/05MR0TEm7sspJlVRLuODbUM9d6ZTaMb44TRZFmz0yStlfQvyYWg\nU3kb8DvAX0fEPZIOBL6UpUAde9FExH/QaMxfCiBpL+C3gN0l7R4R92V5kV5wTxqz/hrGkO90odOT\nm+7k13fe1XHzFsua97gK+HJEbJP0ThpX8h/faacR8VPgXan5e4CPddpmQqZuksm4B5+kccXpw8A8\nYAONdiYzGwJDcSeyDv3gZ84/mJnzD56cf2z1N5tXGQXmpuZ3usan6Qbbn6MxLldHkg4GPgq8mEYF\ne2JfL5hq26wnWT8CHAXcEREH0jji3JhxWzOroKH8xZyvm+QIMF/SPEm7AUto1NgnSdonNXsqOw7C\n2M7FwP8FngKOAy4lYxNN1oDfFhG/AGYkV6BeB7w847Y9U/vag1mJDOP3LU8vmogYA5YBa4D1NHrL\nbJC0QtLERUvvknSbpFuSdd+aoVgzI+LbgJITuB8GXpPl/WS9kvUxSbsD3wUuk/Qw8KuM2/aU2+LN\n+mMov2c5r2SNiNXAgqZly1OPJ6/u78Kvky6YmyQtA34G7J5lw6w1+FOBJ4H3AquBu2j0pjEzq48S\njgcPvAd4Bo0TrUcAb6bRZXJKmQI+In4VEWMR8VREXBIRn06abAZiGH86mvXTUJxQbSFPE02BAvgi\njfb8lwMvpHGCdkodm2gk/ZLWxyoBERHP7q6cZlZmE6E+0TwzdE2iJRpNMuUyGvfUuJUuR8GZqh/8\ns3IUyswqJh3szYaiRl+u0SQnPBIRq6ZebWeDGk0yt6H4sJkNSLrWPkzfNY1nn/pouaTPS1oq6bSJ\nKcuG/bhlX2GG7uejmRWqpLfsextwCI3xZyYOLQF8daoNKx3w4JA3K9pQnXAtZ8D/dkQsmHq1nVW2\niWaCw93Meqac3SRvkPTi6WxY+YAfmpqF2QAMVe2d0naTPApYK+n2ZLj2WyVlGq698k00ZmY1t3i6\nGzrgzWwHw1Rj30kJ2+CTe7FOiwPezHaQ7gs/bGHf5+6Phat8GzwMeY3DrABD+50q50nWaatFwMNw\n1jbMinLN/euGsoda3pOskhZL2ijpDklnt30d6XRJ45IOL+q9QI0CHtxl0sxyylGDT4b0vQBYRONu\nd0slHdJivd2BP6UPN02qVcCbmeWRswZ/JLApuSnHNmAljaHWm/0VjVv1/aawN5JwwJuZTcjXBj8b\n2JKaH02WTZJ0GDAnIq7ucclbcsCbWVvD1uyZc7CxVmMNTx4KJAn4FHDWFNv0TK26SXpcGrPeGrqO\nCx16x/zqvjv51X13dtp6FJibmp8D3J+afxaNtvnvJGG/D3ClpFMi4ubpFrmTWgW8mfXO0IU7dAz4\nZ+4/n2fuP39y/uffX9O8yggwX9I84AFgCbB0ctcRjwN7TcxLug54X0Tc0oOSt1S7Jpqh/FCaFWAY\nfw3nOckaEWPAMmANsB5YGREbJK2Q9NoWLxe4icbMBqXTHZ5qKecFTBGxGljQtGx5m3V/L9+rTa12\nNXgYog+jmfVUSUeTnLZCA17S0yXdJOmWZIjL5cnyiyXdnSy/WdJLe/3aDnmz3hiq71LNhiootIkm\nIn4j6biIeELSLsD3Ja1Onn5/REx5y6k83KvGLJ+hCneqUzPPqvA2+Ih4Inn49OT1JnqQFnpywcym\nb9iCfVLNAr7wNnhJMyTdAjwIfDMiRpKnPiJpraRPSNq16HKYmU2lbm3w/ajBjwMvk/Rs4GvJvQXP\niYiHkmD/HHA28JFW298V6ycf78nzmaW9Wq3WlptpzOppazzMozzS251WJLiz6ls3yYh4XNL1wOKI\n+GSybJuki9nx0t0dHKSF/SqimSWqcC/WWdqLWduvG+Ke2JB/pzUL+KJ70TxP0h7J45nACcBGSfsk\nywS8HrityHKYWXfKHu5FcRNNd/YFLknGSZ4BXB4RV0v6tqTn0TjRuhb4o4LLYWYZDGuwT6pIcGdV\ndDfJW4Gd7lgSEccX+bpm1r2hD3fqd0/WoRiqwCdazdpzsG9XlaaXrIYi4KH1h9ihb8PO4d7EAW9m\nVlM1C/haDjaWlWsvVmX+/PZe3l40khZL2ijpDklnt3j+nZJ+kozD9d1WN+XupaGvwae/JG6ysSrJ\n+3n1AaKFHDX4pLfgBcDxNO7kNCLpyojYmFrtsoj4bLL+62jcwu+k6b9qZ0Mf8GbDxKHemSJXG82R\nwKaI2AwgaSVwKjAZ8BHxn6n1d2f72FyFGOommlb8BbC68md7ajlvuj0b2JKaH02W7fga0h9LuhP4\nGPCuIt7HBAd8ysQXwF8EsyGVbzz4ViPk7rRmRPyfiJhPYwyuP+9FsdtxE00b7jtvVdTuc+tKSzad\n+sE//uCd/PLBuzptPgrMTc3PodEW387lwIVdFK9rDnizGpkYJCwd8g73LnQI+GfvPZ9n7z1/cv7+\ndWuaVxkB5kuaBzwALAGWpleQND8i7kxmXwvckbvMHbiJpgN/MayK/Mtz+vJ0k4yIMWAZsAZYD6yM\niA2SVkh6bbLaMkm3SboZeA9wZpHvxzX4KbipxmyI5LzQKSJWAwuali1PPX5PvlfojgPerMb8K7Q7\nHovGzCrB4d49jdcr4R3wZjXkcJ+meuW7Az4Pj1BpZeRwnz6PBz+E/IWxqvBnNaea1eDdTbLHFu13\nqL9k1nf+zPWG78lqmUx84dxkY/3ikO+BfIONlY4DvmAOerPqcBu8TYuD3orimnvvVKXpJSsHfJ85\n6M1KzE001gtFBP1UNTkfVMw6q1sN3r1oBqyfP6/9U75+/H/aY/nGgy8d1+BLwM02ZuXgGrwVph99\n6F3jM+tgPLJPFeAafAn1I+T9a6H6fLDuvbzdJCUtBs6nUXm+KCLOa3r+vcAfANuAR4C3R8SWnXbU\nI67Bm5lNiMg+NZE0A7gAWAQsBJZKOqRptZuBIyLiMOAK4O+KfDsO+CHl2p/ZznIOVXAksCkiNkfE\nNmAlcGp6hYi4PiJ+nczeCMwu8v044IeYQ96sSb5eNLOBdHPLKJ0D/B3AN/IWuRO3wQ85t8ebbacO\nFzo9uvVuHnv07o6bt1jWcoeS3gwcARzTTfm6VWgNXtLTJd0k6RZJt0paniw/QNKNkm6X9M+SfKAZ\nINfkp89/u5oZbz/t+ZwXcOCBJ0xOLYwCc1Pzc4D7m1eSdAJwLvC6pCmnMIUGfET8BjguIl4GHAac\nJOkVwHnAJyJiAfAYjZ8qNkAOqt7z0NHVo4jMUwsjwHxJ8yTtBiwBVu2wf+llwIXAKRHxi6LfT+E1\n54h4Inn49OT1AjgOWJosvwT4MPDZostillVz01WnoJ5Y12FeAzn6t0fEmKRlwBq2d5PcIGkFMBIR\nVwF/CzwT+FdJAjZHxOt7UPKWCg/4pOvQj4GDgH8A7gIei4iJHqejwH5Fl8Msq4mg7iawpzoAgK9U\nroK8V7JGxGpgQdOy5anHr873Ct0pvBdNRIwnTTRzaHQjelGr1Youh1lWE0F8zf3r+jYYnGv/JZGj\nH3wZ9a2bZEQ8DlwPHAU8J6nZQ5sTERPuivWT09Z4uA8lNduxtl10yE/nFwP4F8HWeHiHfOgFjWef\nqqDoXjTPk7RH8ngmcALwU+A64I3JamcCV7bbx0FaODnN0l5FFteGyFRhWmSNOr3v5tdxTT67Wdpr\nh3zoiZrV4Itug98XuCSprc8ALo+IqyVtAFZK+ivgFuCigsthGXTTJ745iMpYm0yXsVX5Or3frCdY\ne1G2bsplBatGbmdWaMBHxK3A4S2W3wO8osjXtuK0a2YoQyi1C870ic50r5eylHu63Huntzpd6FRF\nHqrAdjDIpovpmgjqLGVr1d491UFhEMr4dx4KY5F9qgAHvO1kuuFSx1CaqOGXrZZfx791GeS80Kl0\nHPDWUlVCvlevl+XEZ79DvtN7K9sBpzZqdpLVAW9tNTd7ZA3TqtYu0++3LL1bsryuw76HahbwHuTL\nLKOqHrisCxXp356VA96mVNZgK2u5eq3qPX2qpCpt61k54K1yhiXYu+Hukj1Ss4B3G7wVotdhM93L\n+c26Mj6efaoAB7xVhsPdCtfhhh87TS1IWixpo6Q7JJ3d4vlXSfqxpG2STivmTWzngDcruSzt7z74\n9UaefvDJkCwXAIuAhcBSSYc0rbaZxvhblxX9XsABbxXg8JqaT8L2SL5ukkcCmyJic3IrvpXAqTvu\nPu6LiNvo06g3PslqhenFBUIOd+urHHd0AmYDW1LzozRCf2Ac8GYV1dx90j1peqBDL5pfPHEfW5+8\nr9PWarXHvEXKwwFvpeOQ2lmrvvBulilAh4B/7sz9ee7M/Sfn79p6Q/Mqo8Dc1HzHmxn1g9vgrXC+\nU1E+vtCpj8bGs087GwHmS5onaTdgCbCqw6u1qvH3lAPe+sK1cquEGM8+NW8aMQYsA9YA64GVEbFB\n0gpJrwWQ9HJJW4DTgQsl3Vrk23ETjVmNuB0+p5xXskbEamBB07Llqcc/AvZv3q4oDngrHQeUDUy+\nXjSl44A3M5vgsWjMrGz8q6dHajYevAPerAI6BbjDvYfGxrJPFeCAt75xEOXT7sbi7h/fQ67Bm02f\nQz6fduHtUO8RB7yZDYoPkAUbj+xTBTjgrXRcG7VBiRjPPFWBA976zrXQYvnvm4Nr8Gb5OYSslNwG\nb9Yb7XqFWGf+mxWoZt0kfSWrDVw6sNz+boMUFbmZdlauwVupuHaajf9OBXETjVmxHF7Z+O9UgJwn\nWSUtlrRR0h2Szm7x/G6SVkraJOkHkua22k+vDGXAb42HB12EaalquaG6Za9auSdCv2rlThto2XOM\nBy9pBnABsAhYCCyVdEjTau8AtkbEwcD5wN8W+XaGMuAf5ZFBF2FaqlpuqG7Zy17udrX4spe7k0GW\nPcYj89TCkcCmiNgcEduAlcCpTeucClySPP4KcHxhb4YhDXgzs5Zy1OCB2cCW1PxosqzlOskdoB6T\nNKuItwLuRWNWG26Tzy/ydX9sdY/V5qp+8zpqsU7PKEp8NlhSeQtnZqUTEdO+kbWke4F5XWzyUETs\nk9r+KODDEbE4mT+nUaQ4L7XON5J1bpK0C/BAROw13TJPpdQ1+Dz/WWZm3YiIA3LuYgSYL2ke8ACw\nBFjatM6/AWcCNwFvBK7N+ZodlTrgzcyqIiLGJC0D1tA4v3lRRGyQtAIYiYirgIuAL0raBPyCxkGg\nMKVuojEzs+mrdS8aSadLuk3SmKTDU8tPkPQjSeskjUg6LvXc4ZJ+klyocP5gSt6x7LMkXSvpl5I+\n3bTNwMtA2KXOAAAEkUlEQVTertzJc+cmF3hskHRiannHi0MGQdJLJd2QfEaulLR76rmW76MMJB2a\nXEBzi6QfSvrt1HOfTsq9VtJhgyxns+Tin5uT6R5JN6eeK+3fu/QiorYTsAA4mEY71+Gp5YcC+ySP\nFwKjqeduAo5MHl8NLCpZ2Z8BHA38IfDppm0GXvYO5X4RcAuNZsEDgDtp9CCYkTyeB+wKrAUOKcFn\n54fAK5PHbwX+Mnn84lbvY9DlTZX7GuDE5PFJwHXJ45OB/5c8fgVw46DL2uE9fBz4UKfPzaDLWJWp\n1jX4iLg9IjbR1DUpItZFxIPJ4/XA0yXtKmkf4FkR8cNk1UuB1/e10NvL2K7sT0TEDcBv0svLUvZ2\n5aZxgcfKiHgqIu4FNtG4MCTLxSGD8MKI+F7y+FvAG5LHp9D6fZTFOLBH8vg5wM+Sx6fQ+EwQETcB\ne0jau//Fy+RNwJeTx+0+N5ZBrQM+C0mnA7ck4TKbxsUJE1pdqFBWZS9780UgP0uWZbk4ZBBuk/S6\n5PGbgDnJ43bvoyzeC3xc0n00LoM/N1le9nIDIOlVwIMRcXeyqBLlLqvK96KR9E0gXROZuHDgzyLi\n36bYdiHwUeDVqW2bFXcRQo6yt9pdi2WFlH2a5W5XvlaVjL6c+e/0PoC3A5+R9BfAKuC/Uus062tP\nhSnKfQLw7oj4elJ5+QKNz3epy5363CwF/rlpnWbuGZJR5QM+Il499Vo7kzQH+CrwluSnHzRqj/un\nVpsD3J+rgB1Mt+xt9K3s0yx3u/IJmNtieeEyvI9FAJIOBl6TLOvrZ6SVTuWW9MWIeHey3lckfT55\nqtTlBkgu/DkNSJ+cH3i5q2yYmmgmawKS9gCuAs6JiBsnlift8o9LOlKSgDOAK/te0p21u+BrcnlJ\ny54u9ypgSTJc6oHAfBonMicvDpG0G41+wav6X9QdSXp+8u8M4EPAhclT7d5HWfxM0jEAko6n0WYN\njXKfkSw/CngsIh4aTBHbejWwISLSAV72v3e5Dfosb5ETjZOMW4AnaVxZ9o1k+Z8BvwRupnGG/mbg\neclzRwC30vhi/H3Zyp48dw/wc+Bx4D6SXidlKPsU5T6XRi+IDSQ9PZLli4Hbk3KfM+jPTVKmdyVl\n2gj8TdNzLd9HGSYaPax+lHyufwC8LPXcBUm515Hq4VSWCbgY+MMWy0v79y775AudzMxqapiaaMzM\nhooD3sysphzwZmY15YA3M6spB7yZWU054M3MasoBb6Uk6ZeDLoNZ1Tngrax8gYZZTg54Kz1Jfyfp\n1uTmG29Klh0j6TpJ/5rcCOKLgy6nWdlUfrAxqzdJbwBeGhEvkbQXMCLp+uTpw2jcgONB4PuSjo7G\nWPlmhmvwVn6/SzJ8bEQ8DHwHmLgN3Q8j4oFojLexlsYdf8ws4YC3smseSTM9n76r1Rj+RWq2Awe8\nldVEkH8X+H1JM5IhfF+Fh4s1y8Q1HiurAIiIryXjl6+jcb/RD0TEw5Je1Gp9M9vOwwWbmdWUm2jM\nzGrKAW9mVlMOeDOzmnLAm5nVlAPezKymHPBmZjXlgDczqykHvJlZTf1/uVbp5jfoExgAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11d847d68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Extract the area variable\n",
    "area = xr.open_dataset('area.nc').load()['cell_area']\n",
    "dom_ds['area'] = area\n",
    "\n",
    "# write the domain file\n",
    "dom_ds.to_netcdf(dom_file)\n",
    "dom_ds.close()\n",
    "\n",
    "# document and plot the domain\n",
    "print(dom_ds)\n",
    "dom_ds.mask.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.close('all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.4"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
